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1.  INTRODUCTION 


The  Rlpt^J  material  response  code  for  x-ray  inter- 
actions includes  three  equations  of  state  for  compacted 
(non-porous)  material:  (1)  modified  PUFF;  (2)  RIP  mixed 

phase;  (3)  GRAY.  A recent  sensitivity  study performed 
by  Systems,  Science  and  Software  (S*)  indicated  that  sub- 
stantial differences  in  radiation-induced  impulse  were 
calculated  depending  on  which  equation  of  state  was  used. 
Furthermore,  John  J.  Ruminer^^J  of  Lawrence  Livermore 
Laboratories  (l’>,  has  determined  that  mixed-phase  effects 
must  be  properly  modeled  in  order  to  obtain  a reasonable 
degree  of  accuracy  when  predicting  the  response  of  porous 
or  tamped  metals.  In  view  of  Ruminer's  results  and  those 
of  another  study  by  Oswald,  McLean,  Schallhorn,  and  Oldham 
at  Harry  Diamond  Laboratories  (HDL) , which  demonstrated  the 
necessity  of  accounting  for  the  melt  regime  when  modeling 
material  response  to  pulsed  energy  deposition,  it  was  de- 
cided that  the  GRAY  three-phase  equation  of  state  was  the 
most  appropriate  for  describing  the  response  of  metals  to  a 
wide  variety  of  radiation  threat  conditions. 

GRAY  is  a three-phase  equation  of  state  for  metals 
developed  by  E.  B.  Roycet^l  at  l’.  Material  near  normal 
density,  in  the  solid-liquid  region,  is  described  by  a 
scaling  law  equation  of  state  for  metals  developed  by 
Grover. [6]  The  scaling  law  equation  of  state  includes  a 
Gruneisen  description  of  the  solid! referenced  to  the 
experimental  Hugoniot.  Material  in  the  liquid-vapor  region 
is  treated  according  to  an  equation  of  stete  developed  by 
Young  and  Alder; t8]  their  van  der  Waals  model  uses  an  analytic 
representation  of  the  classical  hard-sphere  equation  of  state 
with  a van  der  Waals  attractive  term  added.  The  complete 
equation  of  state,  as  described  in  Reference  5,  is  developed 
by  analytically  joining  the  Grover  scaling  law  and  the 
Young-Alder  models  at  a volume  in  the  range  1.3  to  1.5  times 
normal  volume.  This  is  accomplished  by  adding  correction 
terms  to  the  Young-Alder  equation  of  state  emoloyed  on  the 
low  density  side  of  the  join  volume.  Young!’’  recently 
introduced  a modification  of  the  original  GRAY  equation  of 
state  which  permits  calculation  of  the  pressure  in  the 
mixed  phase  liquid-vapor  region;  however,  the  basic  physical 
models  and  the  join  procedure  were  not  altered. 

A substantially  modified  GRAY  equation  of  state  has 
been  developed  and  is  presented  in  this  report.  A critique 
of  the  original  GRAY  equation  of  state  which  illustrates  the 
need  for  improvement  is  presented  in  Section  2.  A continuous 
model  for  the  liquid  phase  is  discussed  in  Section  3.1  and  a 
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soft  sphere  model  for  the  vapor  region  is  developed  in 
Section  3.2.  The  liquid  and  vapor  regions  are  connected 
via  an  improved  join  procedure  which  is  discussed  in 
Section  3.3.  The  results  of  calculations  performed  with 
the  modified  GRAY  equation  of  state  are  presented  in  Sec- 
tion 4 and  a summary  is  given  in  Section  5. 
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2.  CRITIQUE  OF  THE  GRAY  EQUATION  OF  STATE 

The  latest  l’  version  of  GRAY  is  a complete  three- 
phase  equation  of  state  for  metals  and  was,  until  recently, 
the  most  sophisticated  equation  of  state  available  in  the 
RIP  code.  However,  a quick  inspection  of  the  procedure  used 
to  join  the  Grover  liquid  equation  of  state  and  the  Young- 
Alder  modified  van  der  Waals  model,  as  presented  in  Reference 
5,  reveals  the  following  problems: 

1.  The  slope  of  the  equation  of  state  is  dis- 
continuous at  the  join  volume. 

2.  A correction  term,  quadratic  in  temperature, 
added  to  the  caloric  equation  of  state  in  the 
van  der  Waals  region  prevents  the  model  from 
approaching  the  correct  ideal  gas  limit  at 
high  volumes. 

3.  At  high  temperatures  the  correction  terms 
dominate  the  van  der  Waals  model  in  the 
vicinity  of  the  join  volume. 

The  first  two  of  these  problems  are  due  to  the  rather  crude 
nature  of  the  join  procedure.  However,  the  third  problem  is 
more  fund2unental  because  it  indicates  that  the  discrepancy 
between  the  two  models  is  actually  quite  large,  even  in  the 
region  of  the  join  volume. 

The  primary  source  of  the  discrepancy  between  the 
Grover  and  Young-Alder  models  is  the  dramatically  different 
behavior  of  their  atomic  components  (related  to  thermal 
motion  and  ions  of  atoms)  of  pressure  and  energy.  This  dis- 
crepancy is  best  illustrated  by  comparing  the  atomic  compo- 
nents of  the  specific  heats  at  constant  volume  and  the  atomic 
components  of  the  Gruneisen  coefficients.  The  scaled  tem- 
perature dependence  of  the  specific  heat  divided  by  the  gas 
constant  and  the  scaled  temperature  dependence  of  the 
Gruneisen  coefficient,  both  as  represented  by  the  Grover 
liquid  equation  of  state,  are  presented  in  Figure  1 and  Fig- 
ure 2,  respectively.  The  notation  employed  in  these  figures 
is  defined  below: 

T ■ temperature 

Tjjj  ■ melt  temperature  (a  function  of  specific  volume) 

■ specific  heat  at  constant  volume 

R ■ gas  constant 

Y ■ Gruneisen  coefficient 

The  specific  heat  variation  shall  be  discussed  first. 
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Note  from  Figure  1 that  the  specific  heat  in  the 
Grover  liquid  equation  of  state  has  a constant  value  of  3R 
in  the  solid  region,  drops  discontlnuously  to  2.85R  at  the 
melt  temperature,  decreases  linearly  to  3R/2  in  the  liquid 
region,  and  has  a constant  value  of  3R/2  in  the  hot  liquid 
region.  However,  the  Young-Alder  model  has  a temperature 
Independent  specific  heat  of  3R/2  for  all  volumes.  The 
join  volume,  Vj,  employed  in  the  present  version  of  the 
GRAY  equation  of  state  is  the  zero  pressure  intercept  of  the 
boundary  between  the  so-called  liquid  and  hot  liquid  regions. 
This  particular  join  volxime  was  chosen  for  two  reasons: 

(1)  Continuity  betwe<»n  the  Grover  and  Young-Alder  models 
could  be  achieved  by  adding  a single  set  of  correction  terms 
to  the  Young-Alder  model.  If  a smaller  join  volume  were 
selected  it  would  be  necessary  to  add  a second  set  of  cor- 
rection terms  to  the  Young-Alder  model  in  order  to  insure 
continuity  with  the  liquid  region  of  the  Grover  equation  of 
state  at  positive  pressures;  (2)  the  specific  heats  of  both 
the  Grover  and  Young-Alder  models  are  3R/2  at  the  selected 
join  volume.  Unfortunately,  at  the  volume,  Vj,  there  is  a 
large  discrepancy  between  the  Gruneisen  coefficients  of  the 
two  models;  this  point  is  discussed  in  the  following  para- 
graph. 


The  Gruneisen  coefficient,  y,  in  the  Grover  equation 
of  state  is  plotted  as  a function  of  the  normalized  tempera- 
ture, T/Tm,  in  Figure  2.  The  small  variations  of  the  solid 
or  lattice  Gruneisen  coefficient,  Ysr  with  volume  are  not  in- 
cluded in  Figure  2;  additionally,  the  fact  that  the  value  of 
the  solid  Gruneisen  coefficient,  Yor  ambient  volume  and 
temperature  is  material  dependent  (Yq  ~ 2.0  ± 1.0)  is  ignored 
for  convenience  of  presentation.  The  significant  observation 
depicted  in  Figure  2 is  that  y is  a decreasing  function  of 
the  scaled  temperature  until  the  boundary  between  the  liquid 
and  hot  liquid  region  (T/Tm  ® 10)  is  reached;  a slope  dis- 
continuity occurs  at  that  point  and  the  ideal  gas  value, 

Y * 2/3,  is  maintained  for  further  increases  in  scaled  tem- 
perature. The  behavior  of  the  Gruneisen  coefficient  in  the 
Young-Alder  model  is  totally  unrelated  to  that  in  the  Grover 
equation  of  state.  In  particular,  the  Young-Alder  Gruneisen 
parameter  is  independent  of  temperature  and  the  functional 
form  of  the  volume  dependence  is 


(1) 
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where  v is  the  volume  and  Vb  is  the  hard  sphere  volunie.  At 
the  join  volume  currently  used  in  the  GRAY  equation  of  state, 

Y has  a constai.t  value  of  about  4,6  (there  is  some  variation 
with  material)  which  is  more  than  twice  the  maximum  value  of 

Y occurring  in  the  Grover  scaling  law. 

On  the  basis  of  the  critique  of  the  GRAY  equation  of 
state  presented  in  the  paragraphs  above  it  was  decided  that 
several  modifications  would  lead  to  substantial  improvement. 

In  Section  3,1  a new  fom  of  the  Grover  liquid  equation  of 
state  is  developed;  the  new  form  is  mathematically  continuous 
from  T/Tin  = 1 to  T/Tj„  -*•  “ where  the  ideal  gas  limit 
(Cv  ~ 3R/2,  Y = 2/3)  is  reached.  The  advantage  of  the  new 
functional  form  is  that  the  artificial  distinction  between 
the  liquid  and  hot  liquid  regions  is  eliminated  while  pre- 
serving the  Grover  scaling  law  in  the  low  temperature  range 
where  it  is  meaningful  to  within  a few  percent.  Also,  addi- 
tional flexibility  is  gained  with  respect  to  selection  of  the 
join  volume.  In  Section  3,2  incorporation  of  a soft  sphere 
model  into  the  Young-Alder  region  is  discussed  and  the  re- 
quired thermal  and  caloric  state  relations  are  developed;  the 
advantage  of  this  modification  is  that  the  Gruneisen  para- 
meter approaches  the  correct  ideal  gas  limit  at  any  volume 
provided  the  singularity  in  Equation  (1)  is  avoided.  Addi- 
tionally the  location  of  the  liquid-vapor  mixed  phase  boundary 
may  be  correlated  with  experimental  data  through  judicious 
selection  of  the  temperature  dependence  of  the  soft  sphere 
diameter.  A new  join  procedure  which  is  continuous  in  slope, 
thermodynamically  consistent,  and  free  of  incorrect  limiting 
behavior  is  presented  in  Section  3,3. 
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3.  THE  MODIFIED  GRAY  EQUATION  OF  STATE 


3.1  A CONTINUOUS  FORM  OF  THE  GROVER  LIQUID  EQUATION  OF 
STATE 


The  Grover  liquid  metal  equation  of  state  is  developed 
and  discussed  in  Reference  6.  The  equation  of  state  is  based 
on  the  observations  that  a corrected  entropy  of  melting 
(2.32  cal/g-atom  deg)  can  be  considered  constv,  for  all  metals 
and  that  the  liquid  specific  heat  has  a universal  linear 
dependence  on  the  scaled  temperature  (T/Tjjj) ; both  of  these 
observations  are  supported  by  experimental  data.  The  Linde- 
mann  law  which  establishes  the  volume  dependence  of  the  melt- 
ing temperature  is  used  in  coniunction  with  the  above  obser- 
vations to  develop  an  approxit.-»te  equation  of  state  for  liquid 
metals.  However,  as  noted  in  Section  2,  the  linear  scaling 
of  the  specific  heat  breaks  down  at  T/Tm  = 10  where  the  ideal 
gas  limit  is  reached.  For  scaled  temperatures  above  10,  an 
ideal  gas  equation  of  state  with  additive  "effective"  pressure 
and  internal  energy  components  at  zero  temperature  (determined 
by  continuity  requirements)  is  introduced.  In  the  following 
paragraphs  a modified  form  of  the  Grover  liquid  equation  of 
state  is  presented.  The  modification  is  based  on  the  intro- 
duction of  a func  '.onal  form  of  the  liquid  specific  heat  which 
has  a hyperbolic  c-ependence  on  the  scaled  temperature.  This 
new  functional  form  closely  approximates  the  Grover  scaling 
law  in  the  lower  range  of  temperatures  for  which  experimental 
correlation  is  available  and  approaches  the  ideal  gas  limit 
continuously;  thus  the  need  for  a somewhat  artificial  "hot 
liquid"  region  is  eliminated.  The  procedures  employed  to 
develop  the  new  equation  of  state  are  identical  to  those  pre- 
sented by  Grover  in  Reference  6 except  that  the  functional 
form  of  the  liquid  specific  heat  is  hyperbolic  rather  than 
linear. 


The  specific  heat  employed  by  Grover  has  the  form 


= 3R-a(T/Tj^)  , a = 0.15R  , T > T^^^ 


(2) 


Because  Equation  (2)  becomes  unrealistic  at  T/T^^  > 10  we 
select  instead  the  following  relation; 


= 3R/2 


1 + 


(“  i-  * 

' m ' 


-1' 


a = .1,  T > T 


m 


(3) 


Equations  (2)  and  (3)  are  plotted  for  comparison  in  Figure  3. 
Note  that  for  scaled  temperatures  below  4.0,  the  mciximum 
value  for  which  data  are  reported  in  Reference  6,  the  hyper- 
bolic form  departs  from  the  Grover  scaling  law  by  no  more  than 
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6 percent.  At  higher  temperatures  the  hyperbolic  form  has 
the  advantage  of  approaching  the  ideal  gas  limit  continuously. 

In  order  to  obtain  an  expression  for  the  entropy  in 
the  liquid  region  we  employ  the  therroodynairic  relation 


dS 

It 


(4) 


Substituting  Equation  (3)  into  the  above  expression  and  inte- 
grating with  respect  to  temperature  yields 


= 3R  Hn  T - in  |a  |-  + ij  + f (v) 


(5) 


where  is  the  entropy  in  the  liquid  region  and  f (v)  is 
an  undetermined  function  of  volume.  Since  the  specific  heat 
in  the  solid  is  3R,  Equation  (4)  also  yields 


S = 3R  in  T + h(v) 
s 


(6) 


where  Sg  is  the  entropy  of  the  solid  and  h(v)  is  another  un- 
determined function  of  volume.  Solving  Equation  (6)  for 
3R  in  T and  substituting  the  result  into  Equation  (5)  yields 


+ 1 + f(v)  - h(v) 


(7) 


This  expression  for  the  liquid  entropy  is  consistent  with 
experimental  data  on  the  entropy  of  melting  presented  by 
Grover [61  if  we  take 


f(v)  - h(v)  = As 


(8) 


where  As  is  the  "corrected"  entropy  of  melting.  Thus,  the 
entropy  in  the  liquid  region.  Equation  (7) , becomes 


= ®s  - r ^ i)  * 


(9) 
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Following  Grt^er's  procedure  we  obtain  the  entropy 
in  the  two-phase  melt  region,  S^jj,  by  adding  the  solid  and 
liquid  components; 


Sm  = t (l-v)S^ 


(10) 


where  v is  the  liquid  fraction  in  the  melt  region  which  is 
assumed  to  be  linear  in  T. 


V = 


T - T (V) 
s 


AT(v) 


where 


AT(v)  = Tj^(v) 


T^(v) 


is  the  temperature  change  for  melting  at  constant  volume, 
T£(v)  is  the  liquidus  temperature  and  Tg(v)  is  the  solidus 
temperature.  It  should  be  noted  that  Equation  (10)  is 
based  on  the  following  definition  of  the  melt  temperature 


T + T„ 
m _ _S I 

^m  2 


Substitution  of  Equation  (9)  into  Equation  (10)  yields 


m 


> 

+ 

(0 

II 

AS  - Jin  a 1 

i)l 

s 1 

2 I 

(11) 


Noting  that  in  the  melt  region  T A Tm,  we  obtain  the  follow- 
ing approximate  expression  for  the  melt  entropy 


S„  = S + v(AS  - .143R) 
m s 


(12) 


With  Equation  (12)  and  (9)  for  the  entropy  in  the 
melt  and  liquid  regions  respectively,  the  following  thermo- 
dynamic relation  may  be  solved  for  the  Helmholtz  free 
energy.  A: 


. . iih\ 


(»). 


(13) 
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Having  obtained  the  Helmholtz  free  energy,  expressions  for 
the  energy  and  pressure  are  trivial  to  derive.  The  melt 
region  will  be  treated  first. 

Equation  (13)  may  be  integrated  to  obtain 


A = Eq  (v)  ~ J" SdT 


(14) 


where  Eq (v)  is  the  cold  compression  energy.  Substituting 
Equation  (12)  into  the  above  relation  and  performing  the 
indicated  integration  yields 


X 

-/=  gdT  - y v(AS  - .143R)dT 


or,  combining  the  first  two  terms 


% = ^s  - / v(AS  - .143R)dT 
T_ 


(15) 


where  is  the  Helmholtz  free  energy  of  the  solid.  Since 


V = 


T - T 

c 

“St 


(16) 


Equation  (15)  becomes 


143R) 


/ (if 


St 


dT 


which,  after  integration,  yields 


Ajjj  = Ag  - (AS  - .143R) 


T^ 

2St 


TTj 

at" 


S 


J T. 


(17) 
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Expanding  the  expression  in  brackets  and  employing  Equa- 
tion (16)  yields  the  following  expression  for  the  free 
energy  in  the  melt  region: 


\ = As  - (AS  - .143R)  (18) 

The  energy  in  the  melt  region  is  immediately  ob- 
tained from  the  definition  of  the  Helmholtz  free  energy, 

A s E-TS,  and  Equations  (12)  and  (18): 


= 

m 


A + TS  + V 
S s 


.143R) 


or,  more  concisely 


+ 

m s 


.143R) 


(19) 


where  E is  the  extrapolated  energy  in  the  solid. 

The  pressure  in  the  melt  region  is  obtained  from  the 
thermodynamic  relation 


Substitution  of  Equation  (18)  into  the  above  equation  yields 


(4S  - .143R) 


(21) 


or 


"■n,  'fs  * •''S  - •!«'')  37  (^) 


(22) 
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Expressing  v in  terms  of  temperature  and  evaluating  the 
partial  derivative  in  Equation  (22)  yields 


3V  \ 2 / dV  “ 2 dV 


(23) 


Noting  that 


dAT  ^ 


0 


and 


dT 


m 


dV 


Equation  (23)  becomes 


d 


( 


) 


-V 


(24) 


Now,  the  volume  dependence  of  the  melt  temperature  is  defined 
according  to  the  Llndemann  law: 


d in  T 

^ ^ TTnV- 


(25) 


An  alternate  form  of  this  law  is 


X 


in 


(26) 


Substituting  this  definition  into  Equation  (24)  yields 


3 


( 


) 


(27) 
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Substitution  of  Fixation  (27)  into  Equation  (22)  yields 
the  final  expression  for  the  pressure  in  the  melt  region: 


T 

^m  = ^s  " ^ 


(28) 


In  summary,  the  relations  defining  the  new  form  of  the  metal 
equation  of  state  in  the  melt  region  are: 


m 


\ 


m 


m 


+ v(AS  - .143R) 
s 


As  - (AS  - .143R) 


Eg  + V T - {^S  - .143R) 

T 

Ps  - vX  ^ (AS  - .143R) 


(29) 

(30) 

(31) 

(32) 


Following  essentially  the  same  procedure  applied 
above  to  obtain  the  melt  region  equation  of  state  yields 
the  improved  form  of  the  Grover  equation  of  state  in  the 
liquid  region.  Beginning  with  Equation  (9)  for  the  entropy 
in  the  liquid  region  and  applying  the  appropriate  thermo- 
dynamic relations  in  sequence  yields: 


= Sg 


- p in  (a  |-  + ij  + AS 


(33) 


3RT 


m 


2a 


AT 

5" 

(AS  - 

.143R)  - (T 

- ■'n,>  (‘S  * 

(“ 

in  [a  |-  + 

1)  - (a+l)An(a+l) 

(34) 
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E.  = E„  + T LS  + .143R 
x>  s in 


(?) 


3RT. 


m 


2a 


An  (a  |-  + 1 

) - (a+1)  An  (a+1)  - a/| : 

' m 

' \ ■‘’m 

(35) 


- X ^ AS 


3R 

5^ 


X ~ 
^ V 

In 

) - (a+1)  An  (a+1)  - a(|-  - ] 

^)l 

. 

' m 

' ' m 

(36) 


As  mentipned  earlier  in  this  section,  the  above  set  of 
equations  replace  the  two  sets  employed  by  Grover  to  repre- 
sent the  liquid  and  hot  liquid  regions  in  Reference  6.  This 
particular  feature  of  the  improved  form  of  the  Grover  equa- 
tion of  state  results  in  substantial  flexibility  with  respect 
to  selection  of  an  appropriate  join  procedure.  The  intro- 
duction of  a 3oft  sphere  model  into  the  modified  van  der  Waals 
equation  of  s:ate  is  discusi^ed  in  the  next  section. 


3.2  A SOFT  SPHERE  MODEL  FOR  THE  VAN  DER  WAALS  EQUATION 
0^  STATE 

The  classical  van  der  Walls  equation  of  state  has  the 
following  functional  form: 


RT 


V - V. 


a 

"7 


(37) 


where  the  subscript  v Indicates  the  vapor  region.  The  second 
term  on  the  right  of  the  equal  sign  represents  the  effect  of 
a mean  field  attractive  potential  on  the  pressure  in  the 
medium;  a is  the  attractive  potential  coefficient.  The  first 
term  on  the  right  of  the  equal  sign  is  identical  to  the  ideal 
gas  law  except  for  the  appearance  of  V],,  the  hard  sphere  ex- 
cluded volume,  in  the  denominator;  this  term  has  the  effect 
of  Increasing  the  pressure  at  small  volumes  thereby  account- 
ing for  the  repulsive  potential  in  an  approximate  manner. 

The  opposing  forces  Included  in  the  van  der  Waals  model  cause 
the  appearance  of  loops  on  Isotherms  plotted  in  Gibbs  free 
energy  versus  volume  space.  A procedure  for  defining  the 
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liquid-vapor  phase  boundary,  often  referred  to  as  a Maxwell 
constiniction,  may  be  applied  to  the  Gibbs  free  energy.  How- 
ever, experimental  evidence  indicates  that  the  liquid-vapor 
phase  boundary  derived  from  the  van  der  Waals  model  is 
highly  inaccurate. 

A modified  form  of  the  van  der  Waals  model  was  pre- 
sented by  Young  and  Alder  in  Reference  8;  this  modified  form 
was  found  to  be  relatively  accurate  on  the  basis  of  the  cor- 
relation between  theoretical  and  experimental  critical  points 
of  metals  and  was  therefore  selected  to  represent  the  vapor 
region  in  the  GRAY  equation  of  state.  The  new  equation  of 
state,  also  referred  to  as  the  Young-Alder  model,  is 


P = — 


RT  1 + n + a 

5 ? 

(1  - n)'^ 


n = ^ 


(38) 


(39) 


The  first  term  in  Equation  (38)  is  an  algebraic  approxima- 
tion tiOl  to  the  hard  sphere  pressure  calculated  by  computer 
experiments.  liU  An  ideal  gas  caloric  equation  of  state  with 
a potential  energy  term. 


E = 


(40) 


is  used  with  the  Young-Alder  model  to  complete  the  descrip- 
tion of  the  vapor  region  in  GRAY. 


It  was  noted  in  Section  1 that  the  Gruneisen  coeffici- 
ent for  the  Young-Alder  region  in  GRAY  is  given  by  Equation  (1) 
which  is  repeated  below  for  convenience: 


f ** 


At  the  join  volume  between  the  Grover  and  Young-Alder  models 
in  the  GRAY  equation  of  state , Equation  (41)  yields  an  un- 
reasonably high  value  of  about  4.6  (there  is  some  variation 
with  material).  Additionally,  the  Young-Alder  Gruneisen 
coefficient  is  independent  of  temperature.  These  unrealistic 
characteristics  may  be  eliminated  by  allowing  the  hard  sphere 
excluded  volume,  Vjj,  to  vary  as  a function  of  temperature  such 
that  Vb  ■+•  0 as  the  temperature  approaches  infinity.  Note  by 
inspection  of  Equation  (41)  that  setting  Vb  equal  to  zero 
yields  the  correct  ideal  gas  high  temperature  limit  of 
Y = 2/3. 

Several  techniques  for  relating  the  temperature  depend- 
ence of  the  sphere  diameter  to  the  behavior  of  a fluid  with 
realistic  repulsive  forces  are  available  in  the  open  litera- 
ture; the  selection  of  the  most  appropriate  technique  may 
ultimately  depend  on  the  results  of  experimental  definition 
of  metal  liquid-vapor  phase  boundaries.  A simplified  tech- 
nique for  obtaining  the  desired  relation  is  discussed  in  a 
paper  by  Anderson,  Weeks  and  Chandler. tl2]  The  technique  is 
called  the  generalized  Rowlinson  method  and  provides  that  the 
hard  sphere  diameter,  d,  be  calculated  according  to 


00 

■/ 


U - 


^-u(r)/kTj 


dr 


(42) 


where  r is  the  distance  from  the  particle  center,  u(r)  is  the 
repulsive  potential,  and  k is  the  Boltzmann  constant.  Grover t6] 
has  shown  that  his  liquid  scaling  law  is  in  close  agreement 
with  computer  experiments  on  the  liquid  phase  for  particles 
having  the  inverse  twelve  potential  which  may  be  expressed 
algebraically  as 


u(r) 


^(f) 


12 


(43) 


It  is  therefore  reasonable  to  expect  that  this  potential 
function  is  most  likely  to  minimize  the  differences  between 
the  Grover  and  Young-Alder  equations  of  state. 

Substitution  of  the  inverse  twelve  potential  into 
Equation  (42) , after  some  manipulation,  yields 


d 


(44) 
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for  the  soft  sphere  diameter.  The  excluded  volume,  Vb,  may 
be  obtained  from  the  soft  sphere  diameter  via  the  relation 


" ^av 
6A 


which,  when  combined  with  Equation  (44)  yields 


(45) 


''b  " ;i7T 


(46) 


where 


0556  ea 


12 


1/4 


In  order  to  avoid  the  singularity  at  absolute  zero  in  Equa- 
tion (46)  and  provide  reasonable  values  for  Vb  at  temperatures 
well  below  critical,  we  introduce  the  following  relation; 


''b  = ''b.  - =2^  - <^3'^  ! T i 2T. 


— m. 


(47) 


where 


1.1 


9 C, 


■2  ^ \ 


KT 


5 C 


This  relation  is  continuous  with  Equation  (46)  through  the 
first  derivative  at  twice  Triq,  the  melt  temperature  at 
ambient  volume. 

Before  leaving  this  section , the  complete  "soft 
sphere"  equation  of  state  will  be  developed  on  the  basis 
of  the  following  relations: 


''b  “ 

^ Vj^(T) 

(48) 

Lim 

V.  = 0 

(49) 

D 

Lim 

(ar)  ■ » • 

(50) 

The  thermal  equation  of  state  is  obtained  by  simply  substi- 
tuting Equation  (48)  into  Equations  (38)  and  (39)  yielding 


p * liiL 

* « « 


RT 

V 


- _ a 

(1  - n)^  ^ 


1 + n -*■  n 


n 


Vj^(T) 

V 


(51) 


The  caloric  equation  of  state  is  obtained  from  the  equation 
of  thermodynamic  consistency. 


(52) 


Substituting  Equation  (51)  into  Equation  (52)  yields 


(53) 
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where 


f « f(T,V) 


1 + 


^ 2 

n + ri 

J 

(1  “ n)"* 


(54) 


Equation  (53)  and  the  thermodynamic  relations 


A = E - TS 


(55) 

(56) 

(57) 


may  be  solved  for  the  thermal  energy  E and  the  Helmholtz 
free  energy  A.  The  solutions  are 


^ _ 3RT  2RT^  (2  - ri) 

- — - — ar- 


(58) 


A 


V 


7 - Sn 
. 2(1  - n)''- 


£n  V 


RT 


(59) 


The  sublimation  energy#  Eg#  is  included  in  the  caloric  state 
relation  in  order  to  insure  that  the  zero  energy  reference 
is  consistent  «ith  the  solid,  melt  and  liquid  regions.  The 
energy  is  referenced  to  zero  at  300®K  and  ambient  volume. 

7in  improved  procedure  for  joining  the  modified  forms 
of  the  Grover  and  Young-Alder  models  is  presented  in  the 
next  subsection. 
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3.3 


AN  IMPROVED  JOIN  PROCEDURE 


As  pointed  out  in  Section  2 of  this  report,  the  join 
procedure  in  the  current  version  of  the  GRAY  equation  of 
state  has  a slope  discontinuity  at  the  join  volume,  Vj,  and 
employs  correction  terms  which  prevent  the  Young-Alder 
model  from  approaching  the  correct  ideal  gas  limit  at  high 
temperature  and  volume.  These  problems  may  be  eliminated 
by  applying  a mixing  function  to  the  improved  forms  of  the 
Grover  and  Young-Alder  models  presented  in  Section  3.1  and 
Section  3.2,  respectively. 

In  order  to  assure  thermodynamic  consistency  the  mix- 
ing function  is  applied  initially  to  the  Helmholtz  free 
energy  of  the  two  regions  to  be  joined.  Letting  = fm(V) 
represent  the  mixing  function,  the  Helmholtz  free  energy  in 
the  join  region,  Aj,  may  be  represented  as 

Aj(T,V)  = fj^(V)A^(T,V)  + [l-fjj^(V)]Aj^(T,V)  (60) 


where  £ V £ V2  and  the  quantities  and  Av  represent  the 
Helmholtz  free  energy  in  the  liquid  and  vapor  regions,  re- 
spectively. A schematic  illustration  of  the  various  regions 
in  the  improved  GRAY  equation  of  state  is  presented  in  Fig- 
ure 4. 

The  algebraic  form  of  the  mixing  function  is  somewhat 
arbitrary  except  for  the  requirements  of  continuity  which 
may  be  expressed  mathematically  as  follows: 


(61) 

^m<^2)  = 1 

(62) 

/ df  \ 

V / V=V^  " 

0 

(63) 

/ df  \ 

\ v=v^  “ 

0 . 

(64) 
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The  following  trigonometric  relation  meets  all  of  the  above 
continuity  requirements  and  is  therefore  an  appropriate 
mixing  function. 


1 

2 


I cos 


(65) 


Having  completed  a join  between  the  improved  Grover 
and  Young-Alder  models  it  is  necessary  to  obtain  thermal  and 
caloric  state  relations  which  are  applicable  within  the  join 
region.  This  task  is  easily  accomplished  in  a thermodynami- 
cally consistent  manner  because  the  mixing  function  has  al- 
ready been  applied  to  the  Helmholtz  free  energy  from  which 
all  the  thermodynamic  quantities  of  interest  may  be  obtained. 
The  pressure  in  the  join  region,  P j , is  obtained  from  the 
relation 


(66) 


Substitution  of  Equation  (60)  into  Equation  (66)  yields 

df 

Pj(T,V)  = [A^(T,V)  - A^(T,V)]  ^ + f^P^(T,V) 

+ (l-fj^)Py  (T,V)  . 


(67) 


Similarly,  the  thermodynamic  relations 

S = - (H)^  '«) 

and 

E = A + TS  (69) 


yield  the  following  thermodynamic  functions  for  the  entropy, 
Sj,  and  the  energy,  Ej,  in  the  join  region  when  combined 
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with  Equation  (60) . 

Sj(T,V)  = fS^(T,V)  + (l-f)S^(T,V)  (70) 

Ej(T,V)  = fE^(T,V)  + (l-f)Ej^(T,V)  . (71) 


Thus  the  equation  of  state  in  the  join  region  is  completely 
defined  in  terms  of  the  mixing  function  and  thermodynamic 
state  variables  in  the  liquid  and  vapor  regions. 
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4.  DISCUSSION  OF  RESULTS 

The  modified  GRAY  equation  of  state  requires  several 
pareuneters  not  specified  for  the  version  developed  by  Royce.f^J 
A complete  description  of  the  input  data  generation  procedures 
for  modified  GRAY  is  presented  in  Appendix  A of  this  report. 

Critical  point  constants  for  three  metals  calculated 
from  the  modified  G^Y  e^^at ion  of  state  are  compared  to  esti- 
mates by  Grosse^^^»^^»^^i  and  to  values  calculated  from  the  L* 
version  of  GRAY  in  Table  1,  The  estimates  by  Grosse  were  ob- 
tained in  several  ways.  One  way  is  to  get  T,,,  the  critical 
temperature,  by  matching  the  entropy  of  vaporization  at  some 
Tq  to  a curve  of  versus  T/Tc  obtained  from  mercury  data. 

Then  one  reads  off  the  corresponding  value  of  Tq/T<,  and  cal- 
culates Tc.  Grosse  has  also  measured  the  densities  of  liquid 
metals  from  their  melting  points  to  their  boiling  points. 

These  data,  combined  with  the  law  of  rectilinear  diameter  of 
Cailletet  and  Mathias,  and  T^  obtained  above,  give  the  critical 
density  Pc.  Extrapolation  of  vapor-pressure  curves  then  yields 
Pc.  Alternatively,  T^  can  be  calculated  by  requiring  the  boil- 
ing point  density  divided  by  to  equal  4.35,  and  using  the 
law  of  rectilinear  diameter.  Grosse 's  estimates  of  the  alkali- 
metal  critical  constants  compare  favorably  with  later  experi- 
mental data.  Note  that BRL -GRAY  is  consistently  closer  to  the 
estimated  values  than  the  L*  version.  Since  the  critical 
point  lies  outside  the  join  region,  which  is  confined  to  a 
small  range  of  volume  in  the  modified  GRAY  equation  of  state, 
the  critical  point  constants  are  determined  directly  by  the 
parameters  of  the  soft  sphere  model.  This  characteristic 
greatly  facilitates  the  process  of  fitting  available  critical 
point  data. 

A Maxwell  construction  was  performed  on  aluminum,  iron 
and  uranium  to  obtain  the  thermodynamic  state  variables  at  a 
series  of  temperatures  on  the  liquid-vapor  phase  boundary. 

A Maxwell  construction  involves  calculating  the  Gibbs'  free 
energy  along  a series  of  isotherms.  Points  on  an  isotherm 
having  equal  values  of  Gibbs'  free  energy  are  indicative  of 
liquid-vapor  phase  equilibrium  and  are  therefore  located  on 
the  phase  boundary.  The  data  are  reported  in  Tables  2,  3, 
and  4 for  aluminum,  iron,  and  uranium  respectively.  The 
variable  names  appearing  in  these  tables  are  defined  below: 

TMP  ■ temperature  on  the  line  (*K) 

VMN  **  volume  on  liquid  end  of  tie  line  (cm^/gm) 

VMX  ■ volume  on  vapor  end  of  tie  line  (cm^/gm) 
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EVMN  » energy  at  VMN  and  TMP  (ergs/giti) 

EVMX  = energy  at  VMX  and  TMP  (ergs/gm) 

2 

PRES  s=  pressure  on  tie  line  (dynes/cm  ) 

The  volumes  which  bound  the  join  region  of  BRL-GRAY  are  given 
in  Table  5.  The  procedure  for  reading  these  data  into  modi- 
fied GRAY  is  discussed  in  Appendix  A. 


TABLE  5.  VOLUMES  AT  BOUNDARIES  OF  THE  JOIN  REGION  IN 
BRL-GRAY 


5.  SUMMARY 


Several  improvements  have  been  incorporated  into  a 
modified  version  of  the  GRAY  equation  of  state.  The  liquid 
region  now  exhibits  a continuous  variation  of  the  specific 
heat  from  an  appropriate  value  at  the  melt  temperature  to 
the  ideal  gas  value  of  3/2  R as  the  temperature  approaches 
infinity.  A soft  sphere  model  has  been  introduced  in  the 
vapor  region  yielding  the  correct  ideal  gas  limit  at  high 
temperature.  A new  join  procedure  which  affects  a limited 
volume  range  is  employed  between  the  liquid  and  vapor 
regions.  Critical  point  constants  generated  from  the  new 
equation  of  state  are  in  closer  agreement  with  estimates  by 
Grosse  than  those  generated  with  the  original  version. 
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APPENDIX  A 


DATA  GENERATION  FOR  MODIFIED  GRAY 


Prior  to  making  a hydrodynamic  calculation  with  the 
RIP  code  using  the  GRAY  equation  of  state,  E0S3,  preliminary 
calculations  to  determine  the  specific  volume  limits  of  the 
join  region  and  the  state  variables  on  the  mixed  phase 
boundary  must  be  made  for  each  material  in  the  mesh.  Fig- 
ure A,1  shows  a P-V  plot  of  the  regions  defined  by  the  GRAY 
equation  of  state. 

The  first  set  of  preliminary  calculations  involves 
making  "driver"  runs  to  determine  the  specific  volume  limits 
of  the  join  region,  Vi  and  V2*  As  previously  discussed  in 
the  text  of  this  report,  the  join  region  connects  the  liquid 
and  vapor  phases.  Equations  from  both  of  these  phases  are 
used  in  the  join  region  with  a weighting  function  to  bias  the 
pressure  to  the  most  dominant  phase.  The  weighting  function 
is  trigonometric  in  form  and  the  pressure  equation  in  the 
join  region  is  written  to  maintain  thexrmodynamic  consistency. 
Figure  A. 2 shows  a schematic  diagreun  of  representative  iso- 
therms in  the  join  region.  The  values  of  and  V2  must  be 
chosen  such  that  no  pressure  maximum  and  minimum  occur  along 
an  Isotherm  within  the  join  region.  An  example  of  this  be- 
havior is  isotherm  T3  in  Figure  A. 2.  Maximum  and  minimum 
pressures  are  not  allowed  in  the  join  region  because  errors 
will  be  introduced  in  the  calculations  of  the  mixed  phase 
boundary  state  points  which  follow.  Therefore,  the  deter- 
mination of  Vi  and  V2  becomes  a trial  and  error  procedure 
in  which  an  isotherm  driver  program  is  run  with  trial 
values  of  Vi  and  V2  until  no  Isotherms  are  produced  which 
have  a maximum  and  minimum  characteristic. 

The  Isotherm  driver  may  be  incorporated  temporarily 
into  subroutine  ENCALC  of  the  RIP  code  or  it  may  be  written 
Independent  of  RIP.  The  driver  should  essentially  consist 
of  a double  loop  procedure  in  which  the  inner  loop  increments 
volume  and  the  outer  loop  increments  user  selected  tempera- 
tures. Within  the  inner  loop,  subroutine  ET  should  be 
called  to  calculate  the  appropriate  specific  internal  energy 
for  the  driver  values  of  specific  volume  and  temperature  and 
then  subroutine  E0S3  should  be  called  to  calculate  the  pres- 
sure and  sound  speed  for  the  given  values  of  density  and 
specific  energy.  One  coding  change  is  required  in  sub- 
routine E0S3  when  it  is  used  in  the  isotherm  driver  runs. 

The  mixed  phase  flag,  MIXFLG,  must  be  set  to  0 to  bypass  the 
mixed  phase  region.  A comment  card  exists  in  the  subroutine 
to  indicate  where  this  change  must  be  made.  The  values  of 
Vi  and  V2  are  input  to  the  program  through  the  equation  of 
state  constants  array  in  the  locations  ESTCON  (14, N)  and 
ESTCON  (15,N)  respectively. 
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Figure  A.l.  Schematic  diagram  of  regions  in  the  modified 
GRAY  equation  of  state. 


L 


The  selected  isotherms  should  range  from  approximately 
twice  the  critical  temperature  down  to  300*K.  Driven  specific 
volumes  should  encompass  the  trial  values  of  and  V2.  In 
general,  the  lower  temperature  isotherms  will  display  the 
undesired  pressure  reversal  characteristics.  However,  the 
lower  temperature  isotherms  also  will  produce  negative  pres- 
sures which  will  not  be  important  with  respect  to  the  genera- 
tion of  mixed  phase  boundary  values.  The  user  should  attempt 
to  obtain  isotherms  which  give  positive  pressures  in  the  join 
region  with  no  pressure  minimums  or  maximums  for  as  low  a 
temperature  as  possible.  In  general,  this  minimum  temperature 
will  be  in  the  range  of  1000  to  3000  *K.  Initial  values  of 
Vi  and  V2  should  be  chosen  as  follows.  Vi  must  be  larger  than 
the  specific  volume  at  which  the  llquidus  intersects  the  zero 
pressure  line,  Vr , and  V2  may  be  estimated  by  1.3  x V^.  Sub- 
sequent values  ot  Vi  and  V2  are  based  upon  the  results  of  the 
isotherm  driver  for  the  previous  values  of  and  V2* 

Once  the  limiting  values  of  the  join  region  have  been 
determined  the  remaining  step  is  to  generate  the  mixed  phase 
boundary  values.  These  values  are  calculated  by  the  DATGEN 
and  ET  subroutines.  To  activate  the  DATGEN  subroutine,  a 
change  to  subroutine  ENCALC  is  required.  The  cells  to  DATGEN 
and  EXIT,  which  exist  as  comment  cards  in  ENCALC,  must  be 
activated  by  removing  the  "C"  in  column  1.  In  addition,  the 
mixed  phase  flag,  MIXFLG,  must  be  set  to  0 in  subroutine 
EOS3. 

The  DATGEN  subroutine  is  based  upon  the  "Main  Subrou- 
tine" described  in  Reference  9.  The  mixed  phase  boundary 
state  points  and  the  critical  point  are  calculated  by  this 
subroutine.  Output  is  in  the  following  form.  For  each 
value  of  temperature,  which  is  automatically  selected  in  the 
DATGEN  subroutine,  values  of  specific  volume,  specific  energy 
and  pressure  are  calculated  on  the  mixed  phase  boundary.  Two 
values  of  specific  volume,  and  and  specific  energy, 

Emax  %IN»  produced  at  each  temperature.  The  maxi- 
mum values  correspond  to  the  mixed  phase  boundary  values  at 
the  specific  volume  greater  than  the  critical  point  and  the 
minimum  values  correspond  to  the  mixed  phase  boundary  values 
at  the  specific  volume  less  then  the  critical  point.  Mixed 
phase  boundary  values  calculated  by  subroutine  DATGEN  are 
used  in  the  /GRAY/COMNON  arrays  in  the  hydrodynamics. 

The  equation  of  state  constants  and  mixed  phase  bound- 
ary data  are  loaded  into  the  ESTCON  array  of  /EQST/COMMON 
and  the  /GRAY/COMMON  arrays  respectively.  Table  A.l  lists 
the  required  ESTCON  variables. 
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TABLE  A.l.  THE  ESTCON  ARRAY 


Location 

Symbol 

Definition 

ESTCON (1,N) 

^0 

3 

Ambient  density  (g/cm  ) . 

ESTCON (Z,N) 
to 

ESTCON ( 5, N) 

* 

Unused. 

ESTCON ( 6, N) 

=0 

Bulk  sound  speed  (cm/sec) . 

ESTCON ( 7, N) 

AMU 

2 

Shear  modulus  (dyne/cm  ) . 

ESTCON ( 8, N) 

’'o 

2 

Yield  strength  (dyne/cm  ) . 

ESTCON  ( 9, N) 

0 

2 

Spall  strength  (dyne/cm  ) , 

ESTCON (10, N) 

=S 

Sublimation  energy  (erg/g) . 

ESTCON (11, N) 

=LV 

Energy  at  which  the  matrrial 
commences  to  vaporize  (erg/g) . 

ESTCON (12, N) 

®LM 

Energy  at  which  the  material 
has  completed  melting  (erg/g) . 

ESTCON (13, N) 

®SM 

Energy  at  which  the  material 
commences  to  melt  (erg/g) . 

ESTCON (14, N) 

VI 

Specific  volume  at  boundary 
between  the  join  region  and 
the  liquid  region  (cm3/g) , 

ESTCON (15, N) 

V2 

Specific  volume  at  boundary 
between  the  join  region  and 
the  vapor  region  (cm^/g) , 

ESTCON (16, N) 

S 

Hugoniot  parameter 

U **  C + SU_  (cm/ysec)  . 
s p 

ESTCON (17, N) 

^0 

Lattice  gamma 
X(V)  B Xq  - ax. 

;ON(18,N) 

a 

See  ESTCON (17, N) . 

ESTCON (19, N) 

\ 

e 

Electronic  gamma. 

ESTCON ( 2 0,N) 

Electronic  energy  coefficient 
(mbar  cm^/mole  deg‘) . 

ESTCON (21, N) 

Melting  temperature  parameter 
(deg) . 
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TABLE  A.l  (Continued) 


Location 

Symbol 

Definition 

ESTCON(22,N) 

®OH 

Energy  at  V = Vq,  T = 300®K, 
P » 0 (mbar  cra2/g) , 

ESTCON(23,N) 

E 

00 

Energy  difference  between  Eqh 
and  energy  at  V = Vq»  T = 0 
(mbar  cm^/g) . if  undefined 
E«  = -300  (3X8.134E-5  + 

150  ge)/AW. 

ESTC0N(24,N) 

AW 

Atomic  weight  (g/mole) . 

ESTCON(25,N) 

Vj 

Unused. 

ESTC0N(26,N) 

^b 

Excluded  volume  for  vapor  phase 
(cm3/g) . 

ESTC0N(27,N) 

a 

y 

Coefficient  of  attractive 
potential  for  vapor  [mbar  (cm^/ 
mole) 2] . 

ESTC0N(28,N) 

^crit 

Critical  temperature  (*K) . 

ESTC0N(29,N) 

AS 

Entropy  of  melting  (mbar  cm^/ 
mole  deg) . 

ESTCON(30,N) 

EOSFLG 

Equation  of  state  flag.  Must 
have  the  value  of  3.0  to  indi- 
cate E0S3  is  to  be  used. 

Table  A. 2 lists  the  /GRAY/COMMON  variables  that  are 
to  be  loaded  with  the  mixed  phase  boundary  data  generated  by 
subroutine  DATGEN. 


TABLE  A. 2.  THE  GRAY  COMMON  ARRAYS 


Location 

Symbol 

Definition 

TMP(I,N) 

1-1,20 

T 

Temperature  (*K) . Load  in 
descending  order  starting 
with  the  critical  temperature 

VMX(I,N) 

1-1,20 

V 

MAX 

Specific  volumes  which  are 
larger  than  the  critical  value 
(cm3/g) . Load  in  descending 
order  starting  with  the  critical 
specific  volume. 
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TABLE  A. 2 (Continued) 


Location 

VMN(I,N) 

1=1,20 


EVMX(I,N) 

1=1,20 


EVMN(I,N) 

1=1,20 


PRBS(I,N) 

1-1,20 


The 

in  exactly 
array. 


Symbol 

V 

MIN 


P 


Definition 

Specific  volumes  which  are 
smaller  than  the  critical 
value  (cm3/g) , Load  in  de- 
scending order  starting  with 
the  critical  specific  volume. 

Specific  energies  at  correspond- 
ing specific  volume,  (mbar/ 

cm2) . Note:  DATGEN  ui^s  which 

are  dyne/cm^  must  be  converted 
to  mbar/cm^. 

Specific  energies  at  correspond- 
ing specific  volumes,  (mbar/ 

cm2).  Note:  DATGEN  units  which 

are  dyne/cm^  must  be  converted 
to  mbar/cm2. 

Pressure  (mbar/cm^) . Load  in 
descending  order  starting  with 
the  critical  pressure.  Note: 
DATGEN  units  which  are  dyne/ 
cm2  roust  be  converted  to  mbar/ 
cm2. 


index  "N"  in  the  /GRAY/COMMON  variables  is  treated 
the  same  manner  as  the  "N"  index  in  the  ESTCON 
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APPENDIX  B 


LISTING  OP  xMODIPIED  GRAY 


C 


c 

c 


c 

c 

c 

c 

c 


SUBROUTINE  FOS3 

CSC 

10 

CSC 

20 

COMMON 

AMASS! lOi .OELTMI ( 1 C 1 tOELTNN 1 1 3 1 tNATSMO 1 1 C 1 1 MATL 1 0 1 10 1 i 

CSC 

30 

• 

RATIO! lOi iNHTKLS 

CSC 

SO 

COMMON 

VISCOL!  lot  iVlSCOgllO)*COURST,l>CLTIH.MAXCri.iMOOEiNRCSRT.TNAX 

ESC 

SC 

common 

IPLATE.UZERU 

ESC 

60 

COMMON 

ALPHA! Al ) lAMU! 1 1 | tBBTCMP ! 2 1 • CON VRG .CHCRfir 1 2 I |M AXPAS tNALPH A • 

ESC 

70 

• 

NFLUO! 13) iNPC! ICt .NSCTT! 1C) .SSTARTI2I .SSTOP ! 2 1 iNUNIT 1 . 

ESC 

60 

• 

NUN|T2|NUNIT3«NUnITS.NF)TI 13) iNEDGE«NiNCiNKFILT,NFUNCT 

CSC 

90 

COMMON 

STIMFI  •STIIIC2«ENt<GY,N|0TH 

CSC 

100 

COMMON 

lOUHlTlAl 

CSC 

DC 

COMMON 

NCYEDT!2CI iTED|T!2C) tNABSIS) .mORDISI I iZPlOTlSl iLPLOTIbl I 

ESC 

129 

• 

1 POT, MOMENT  1 NOUMH ,NFEO 1 T «NPK igTiTPRI NT iNSEO IT 

ESC 

I3C 

COMMON 

ICUTLilCUTRilAOOLtlADORtlOOEN.DhZl .0KAT10,A00MAS,0UM| 

ESC 

ISO 

COMMON 

ivzone,viioth 

CSC 

ISC 

COMMON 

OELTaH! IOb?> ,0IScPTII2) ,hERkOR)2) , ITEMPI25) .MZONC! I0,S) , 

ESC 

lt>9 

• 

TCMPI2S) ,uTAVti,DTN,OTNPI .lC0UHN*JSAVEtHAPFL6,hAXZ,MAXZPI , 

ESC 

170 

• 

NCYCLE,NSAVE,PbTH0S,P0TNr6,T|HC 

ESC 

160 

COMMON 

C!  I33;.i  tC!  i:OC)  .HRATEDCSOtZ)  .leFLAfi!  ICCC  1 , 1 VFLAG ! lOCO)  , 

ESC 

190 

% 

OLORHO!  lOCC)  •PlIcCOtOll'CO)  ,RHO  1 1 500 ) ,S  ! 1005 ) ,S0  i I09C  1 • 

CSC 

200 

% 

U!  IC90I tXI IOCS) tULANK 137) 

ESC 

2IC 

CSC 

229 

COMMON 

/COST/  CSTC0N!30,3S) 

CSC 

230 

CSC 

2HC 

COMMON 

/GRAY/  TMPIZCiSOl «VMX!2C*30) .VMN)2C,30) ,EVHX! 20,30) « 

• 

EVMN!20,3C) •PRCSi2C,3C) 

COMMON 

/TEM/  TIIC90) 

COMMON 

/TSAVE/  T0L03,E0t03 

COMMON 

/EC/  EC0N!23,3Ci 

OATA  CAMMAF.ALF/|*C5S6,«|/ 

DATA  B0LTZ«AVN0«P|/|.3«E-U*«*C2E29i9*11|*/ 
C 


C 


C 

C 


ESC 

•••  statement  rUNCTIONS  •••  ESC 

ESC 

EMIXI  •X2iX3l>TENP( 1 1 I♦X|•I3••EC^NIS•(I)«•S•EC0NI2•N)•X| l»X2*  ESC 

IECONI3«NI-*S«CCUNC«l,NI«l  1««X3**2)  I ESC 

CSC 

PC(X4,XSia*S«lCSTCONIIV|NI>CSTCOIll  I 7 |N  I «ESTCONC  iS  .N  I •X'li  •ECOHI  2 ,N  I CSC 
/(CC0NII«NI«IU*X«I))*XS««2  CSC 

CSC 

Tl  (XA.X7iX«l>l-X7«SQRTI«7*«2>2*tXA*XSn/XA  CSC 

CSC 

J-.ISAVC  CSC 

N«NSAVE  CSC 

TOLO>TIJi 

COLO>CIJi 

1J0|N«0 


26C 

270 

260 

270 

3C3 

310 

320 

330 

3««3 

350 

3*0 

370 

380 
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L PRSCSDim  Pi3Ki.BUI0C-l«T  PIIMSI 


IF  leCUNO.NI  tST*  C»)  (iO  TO  3C 

•••  initializatiom  of  ehuation  of  state  constants  ••• 

ECONI I ,N)>I ./ESTCONI I ,N| 

ECON(2tNi>ESTCON(2C«N)/ESTCON(2<**N) 
EC0NI3.NI>CSTC0N(29.N|/t:STC0NI2<t«NI 
EC0N(<t»NI«l.  U3E>G/ESTCUNI2<ltN) 

ECONI5tNI>a.l3‘«E>S/CSTCaN(29tN) 

EC0NI7(N)«‘t.SE9«EC0NI2tNI«ESTC0N( I «N I • (ESTCON 1 1 7 (N l-ESTCON I 19 «N I I 
IF  IESTC0NI23(NI  .6E*  O.i  ESTC0hl23.N)--920*«ll. I39E-S«S0*« 

. ESTCnNI20.rin/ESTC0NI29tNI 

ECONI  |2.NI>l.«ESTCONI  I .Ni«ESTCOM  M«Nl 
ECONI |5(NI«IALF>I.0I«ALUCIALF«I.&I 
TEMPI  I laESTCONI l7tNI*|.i33 
TENPI2IbESTCONI2|  ,NI/(  |.3«CCON(l2,Nn**2 
TEMPI  31 ■! *«FC0NI I2>NI*I2««TEHPM I ♦ECON 1 1 2 .N I • I TEMP  1 1 i« 

. I2*«ESTCUNI |7.NI>|«*«A7l-ESTC0NllfltNI I I 
ECONiatNl«TFMPI2l«TCMP|J| 

TEMP|9la2i«IECONI I2«N|«ESTC0NI IH»NI«ESTCONI I7.NI«*3333I 

TEMPI5I>I»2«I l••CC0NI|2.NI l•l•E•l^•CSTC0NIA,NI••2 

ECONI9.NI>ECUNI3«NI*ITEmPMI«TEMPI2I*TEMP|3M*«2/TEHPISI 

TC>ESTC0NI2A«NI 

VB>ESTCONI2AtNI 

SIGMA>l3*«ESrCUNI29.NI«tfB/(Pl«AVNO*9*n**.33333333 
0T«l  A*  •ESTCONI  29  till*  VD/ 1 PI«AVNO  I I ••*33333333 
EPS>B0LTZ*TC^I0T/I6AMMAF^SI6MAI  MPU 

ECONIAtNl>P|PAVNO«l I GAMmAFpSIGMAp I EPS/BOLTZ I • • *08333333  I pp3  I / 

* IA*9pESTC0NI29*NI I 

ECONI  |3.NI>*l93PECOtllS*NlP*SPECONI9*Nl 
TEMPI  1 laESTCONI U*N)PECUNI I2.NI 

TEMPI2I>*SP| l*E-iPESTCONIA*N|PECONI I2*N) |pp2/I I *«TEMPI I ) > 
TEMPI3)>I*PTEMPI l|P| l*P*SPTEMPI 1 1 • 1 1 *-ESTCON 1 1 7 *N ) /ESTCON I I « *N ) ) I 

* /3» 

ECONI II .N|aTEHFl2IPTEMPl3iPESTCnNI23*N>«l |*PESTCONI |7*N)P 

* CC0NII2*NIIpESTC0NI22*N| 

TEMPI 2Si«TEMPI 21 pCCONIA.NI 

ECONI |7,Ni>ESTCONl27,NI/ESTCONI29*N|P*2 
IF  I ABSITMPI I *NM  *LC«  d*i  GO  TO  22 
00  20  ll>3.20 
JHAX-I !•! 

IF  lABSITMPI I I ,NI I *CC.  0*1  GO  TO  21 

20  CONTINUE 

21  ECONI  |4«NI*VMN  UMAX  *NI«PIIESIJMAX*N  I Pi  VMN  UMAX  >1  ,N  l-VMN  IJNAX  *N  i)  / 

* IPPESUNAX-|,N|«PIIESUMAX*NII 

22  ECONI |9,NI>2*pESTC0NI2I.NI 

ECONI |B,Nla| *|PEC0NIA,N|/IEC0NI|9*N|pp*2Si 

ECONI  |9,N|p2*pEC0NI  |8*Nl/EC0NI|9*NI«2*2SPEC0MU,Nl/ 

* lECONI |9,N|PP|*2SI 

EC0NI20«NI>|*2SpEC0NI4*NI/ICC0NI|<I*NIpp2*2SI-ECUNI lliNI/ 

* lECUNI |9*N|PP2| 

PRINT  2St  lECONI |,N|,|P|,20I 

2S  FORMAT  (//*7H  ECONS •// .9 (9X . IPSE I9.G*/ I . // I 


ESC  390 

ESC  900 
ESC  910 
ESC  920 
ESC  93C 
ESC  99C 
ESC  9SC 

ESC  970 
ESC  990 
ESC  SOO 
ESC  SIO 
EXVI 


TMVI 


TISVI 


TIIVI 


ESC  7S0 


IVMNOI 


IVBOl 

(C2I 

IC3I 


CSC  930 


u u uu 


33  E( J)«I»E-I2*E( J( 

KNU>ESTCON|  ( 

k'X>|.-RNU  / 

•••  FOR  DRIVER  OR  DATGEN  RUNS.  CHANGE  NIXFLG  TO  3*  THIS 
BYPASSES  THE  HIXEb  PHASE  REGION.  ••• 

MIXFI.G>I 

IF  IMIXFLG  .£0.  Cl  GO  Tu  3S 
ITER-0 

V-I.O/RHOUi 
VOS-I.O/ESTCUNM  .til 
IF  IV  .LE.  ECONIU.HII  tiO  TO  3S 
IMX-I 

00  7122  I-I.I1 
K-UI 

IF  (EVHXIK.NI  .6T.  EVHXll.NII  IHX-K 
7122  CONTINUE 

IF  (EIJI  .LT.  I.2«EVHX(IHX*NI(  GO  TO  IIS 
32  EXV2-I.-ESTC0NII .NI-ESTCONIIS.N) 

IF  (EX  .GE.  EXV2(  GO  TO  3S 
TEHPd  l>l.5*£C0N(S.N( 

TENP(2OESTC0N(  I.NI/I  I.C^-EXI 
60  TO  1C 

3S  TEHP(|(«ESTCUNII4.NI«CX 

TEHP(2(«.S«( I.E-4«ESTC0N(4*NI*EX(*«2/( l.-TEHPI I I I 

TEHP(3(-I.*TEMP(I  l•(l•♦•S•TEHP( I )•( I .-ESTCONl I7.NI >ESTCON( I4.NI >1/ 
. 3* 

TEHPd I (>TEHP(2(«TEHP|3|*ESTcONI23.NI*I I .♦ESTCOtd  I 7 ,H I *EX t ♦ 

. ESTC0N(22.N( 

IF  (EX  *GT.  C.l  60  TO  10 
TEHPd  l-ESTCONI  17. Nl-I. 3333 
TEHP(2I>ESTC0NI2I .Nl/I I .-EX(«*2 

TEHP(3(«I«*CX«I2.«TEHP( I l*EX«(TEHPI I l•l^.•eSTCON( I7.NI*I .44i7l- 
. CSTCONdS.NIK 

60  TO  SC 

10  TEHP(|(-2.«(ESTC0(I(I7,Nc.3333> 

TEHP(2(«ESTC0NI2I.NI 

TEHP(3(«I. *EX«( TEHPd |*£X«( (ESTCONI I 7 «NI>.33 33 1 • ( TEHP ( I ) «l . 

• ESTCONdS.MIII 

50  TEHP(|2(«TEHP(2I«TEHPI3| 

TEHPI |1|>2.«(EX«ESTC0II( |i.N(*ESTCON( I7.N)«.3333I 
TEHPd (>|.2«( |.•EX(•I.E-I2•ESTC0N(4.N|••2 
TEHP(2(«|. 

IF  (EX  .6T.  0*1  TEHP(2I«|.»EX*(1.4ESTC0N(I4.NI-|.I 

TEHPI  |SI«EC0N(3.NI«  (TEHPI  l1l*TE«*Pd  21)  ••2/(  TEHPd  l«TEHP  (2d 

TEHPd  («TEHP(  |2I>.S«TCHP(  IS( 

TEHP(|4I«EH(TEHP(|(.(^.,0.I 

TEHP(|(«ESTC0N(|4.NI«CX 

TEHPI2laESTC0N(I.N(/l|.-EXI 

TEHPI3laEX«TEHP(2l«( I.E-4*ESTC0N(4.NI7I I.*TEHP( 1 1 I (••2 
TEHPI1l«l.-EX«(d.«.S-ESTC0Nd7.NI  l-.S*EX*ESTCOM  li.NI  I 
TEHP(|1l>TEHP(3(«TEHP(1|«TeHP(2l*IESTC0N(|7.MI-tR*ESTC0N( li.NI (• 


ESC  110 
ESC  ISO 
ESC  143 
ESC  170 


ESC  1003 
ESCIOlO 
ESCI02C 
ESCIC33 
ESCIOlO 
ESCI0S3 
ESCI040 
ESC«073 
ESCI080 
ESCICIO 
ESCII03 
ESCIIIO 
ESCII20 
ESCII3C 
ESCIIIO 
ESCII50 
ESCI I4C 
ESC  1 1 70 
ESCI I0C 
ESCIIIO 
ESCI20C 
ESCI2I0 
ESCI220 
ESCI230 
ESCI210 
ESCI250 
ESCI240 
ESC  1 270 
ESCI2S0 


S3 


non  n non 


. IC(J)-ESTC0NI22,nI  I 

VB|«/RHO(JI 

IF  IV  .GT.  VII  GO  TO  70 
IF  leiJI  tGr.  TCNFIUII  GO  TO  *7 

•••  SOLID  KCGION  ••• 

IVFLAGI Jl>20 
TEMPI  I laS.BECC^KSiNI 
TEHPI2la-IEIJI"TEHPI I |l  I 

ARGaTEHPI  I l••2•2.•EC0rll2•NI•TEMPI2l 
IF  URG  *LT.  C«l  GO  To  aS 
TEHPISI>T|IECUNI2iNI iTEMPII l•TEHPI2l I 
TIJIaTEMPISI 
TEMPIGIaPCIEXiTEMPISI  I 
PIJIaTEHPI  |9I«TEHPIGUEC0NI7«NI 
GO  TO  IOC 
SS  PlJIaO. 

TIJtaQ. 

GO  TO  100 

60  TEMPI  I taTEHPI |2I«.S«TEHPI ISI 
TEMPI2laTEMPI I I/TEMPI  121 
TEMPI |7|aEMITEMPI| t .TEMPI  121 tTEMP|2l I 
IF  lEIJI  «GE«  TEMPII7II  GO  TO  7? 

•••  melt  region  aaa 
IVFLAGIJIBSI 

TEMPI  I talEI  JI-TEMPI  lAI  l/ITEMPI  I7I«TEHPIUI  I 
TEMPI2laCC0NI3.Nl«EC0IJIt«N| 

TEMPI3la3,aFC0N|S,MI«TEHPI I l•TEHPt2| 

TEMPIRIa-ieUt-TEMPI  1 1 1 ♦•S«TEMPI  t IBaEBTENP  1 1 G I bTEMP  I 2 1 1 
ARG«TEMPI3l**2-2*«ECOi:ii.N|aTEMP|RI 
IF  lARG  .LT.  C*l  GO  TO  GS 
TEMPISIaTI IEC0NI2.Nl .TEMPI3I .TEmPIRM 
TIJiBTEMPISi 

TCMPI2iBESTC0NI|.HI/l|.«EX| 

TCMPI3I«TEMPI5»-.5*TEIIP|  IMTEMPIISI 
TEMPIRIaESTCONI |7.NI-rSTC0NII0.NI«EX 

TEMPI AlaPCIFX.TtMPiSI l-TEMPI I l•TEMP^2l•IEC0Nl3.Ml-CC0NI•I.MI I* 
. ITfMPI  IHtTEMPI  I2I*TEMPI3I*TEHPIR|  | 

PIJiaTCHPI |Ri«TEMPIAI«EC0N|7.NI 
GO  TO  ICO 
C 

7C  TEMPI  I |a|.S*t'CONIS.in 

TEMPI2I>ESTC0NI I .Nl/I I.O-EXI 
TEMPI  A I a. |R3*EC0MIS.N|«.S«TEMPI |S) 

TEMP||CI«2.C*TEMPIII 
TEMPI |3I>TEMPI I I/ALF 
TEMPI22|a,S«EC0NI2.r.'l«reMP|  12) 

TEMP  I 23 1 ■! TEMPI  I I l-E Ul «TEMP I A 1 1 /TEMP  1 1 21 ♦CCUNI 3 .N) ♦ 

. TEMPII3l«IALF-EC0N||GtNI) 


ES02VC 

ESCI3CC 


ESCI3): 

ESCI32S 

ESCI33: 

ESCI3H? 

ESCI3S: 

ESCI3A3 

ESCI370 


ESCI383 

ESCI39C 

ESCMOO 

Escmo 


ESCIR23 

ESCIR33 

ESCIRRO 

ESCIR5C 

ESCIRoO 

ESCM7C 

ESCIRRO 

ESCIRRO 

ESCISOC 

ESCISIO 

ESCIS20 

CSCIS30 

ESCISRO 


CSCISS3 

ESCiStoO 

ESCIS70 

Esciseo 

ESCISV5 

ESCIAOO 

ESCIAIO 

ESCIA20 

ESCIA30 


54 


GAKMASaCSTCOHt  l7tN>-eSTC0NnaiNI*EX 

VI«ESTCON(M«NI 

V2«ESTC0N( IStNI 

I.JOI»l>C 

IF  (V  tLE*  VII  60  TO  7<l 
IF  (V  «6E«  V2I  60  TO  90 
IJOINal 
C 

C •••  JOIN  RE6I0N  **• 

C 

IVFLA6I Jl>33 

TeNF(3l>2.a*EC0N(6*..l*Tt:NP(2l 
TENPMI>-TEHP(2I»EC0N(I7*NI 
TENP(25I>TEHP(2I«EC0«U,NI 
FN>«5*(C0S(P|«(V2-VI/(V2-VI I l♦l•OI 
0FNDVa«S*P|«SIN(Pl*(V2-Vl/(V2-VI I I/IV2-VI I 
C DFNOVbO* 

XV>I lO'ESTCnNI I (NIBV 
TI7aESTC0NII7.NI-l.333 
T2VsESTC0N(2I .NI / ( I .C-XV I«*2 

T3Va|.C«XV*(2.*TI7«XV«(TI7«(2.*ESTC0N( I7.NI-I .*»*7l-ESTC0N( lO.NI 1 1 
THVsT2V*T3V 

T9Va2.«(XV*FSTCON(ia.UI-ESTCON(l7.NI«.3333l 
TSVs|.2aXV«I.C-|2*ESTCON(A*Nla*2 
EC0N9V>EC0N(3.N|a( (T2V«T9V*T3V|a«2l/TSV 
TIVbESTCONI iP.Nlalv 

T2Va.S*( I.E-*«ESTC0N(4*MI*XVI*«2/( I .C-TI VI 

T3Va|,«riVa ( | • « «S*T I V* ( I .-ESTCONI I 7 .NI/ESTCONI U.NI 1 1/3.0 
EOVaT2VaT3V«ESTCON(23.N|a( |•♦ESTCON( |7.N|aXV»«ESTC0N(22.NI 
ITERTjac 

IF  ITOLO  .LT.7.1  TOLOaSoO* 

TJaTOLO 

ITHaO 

TUPaI.EA 

TL«a300. 

73  IF  (TJ  *6T.  ECONIM.NII  60  TO  731 

Ve«EC0N( l•.NI•TJ•(eC0N(|9,NI♦CC0llt2C*NI•TJ| 

TENP(7la-EC0M( I9.NI-2.«EC0N(2C**>I«TJ 
TEHP(6la-2.«EC0NI2S.N| 

60  TO  732 

731  ve«ECON(*.NI/( Tjaa.2SI 
TENP(7la-.2S*ECON(A.N|/|TJ**|.2<il 
TEMP(6|sS.*EC0N(A.N|/( U••TJaa2.2SI 

732  ETAaVB/v 

TENP(Slal2.!}-ETAI/(  ( l.0-ETA|a«3l 

EVaTEHPl I l•TJ-2••CCONIS•NlaTJa•2a{2.-ETAI•TCNP|7l/(V•t I.-ETAI»«3I 
. -ECON(r,7.NI/V«ESTCONllO«N|*|.E-l2 
XXaTJ/TNV 

EL«E6V*TENP( IOI*Tj«.5«CCONI2.NiaTJ*a2«ECONI3.N|aTNV»ECON(9.NI*.$a 
. eC0N9v«TENP( l3l«TNVa(«L06iALFaxx«|.0l-EC0Nt |».Nl-ALF*(XX-l .01  I 
FFa-CI  J|*FH*EV*I  l.-FNiaEl, 

IF  IFF  «6T.  d.l  TUPa«(i|MllTUP.TJI 
IF  IFF  .LT.  0.1  TLN*AHAaI(TLA.TJ| 


V' 


S5 


IF  UTM  .QT.  C>  60  TO  733 

TCMPC II l■TCHP^7»•TEMP^2l  OETAOT 

TENPI  i2|aTCMPU0l*CC0tM^.NI*TJ*TCNP<  1 1 •!  1 1 . /ULFaXXal  •Ol*l  tOOELO? 

TEMP  II II -TEMP U »-1«*ECONl5,N|aTj*TEMP4S»*TEHPI7l/V  (OEVOTI 

• . -2««EC(IN(S.ri)«Tj«TJ*TEMP|SI«TEMP(SI/V 

• ♦2**EC0NIS«N)*Tj*T3«TEHPI7)«TEHP(7I/( V«V*U  t-ETAi^'SI 

• -A««EC0NIS|NI«Tj«T3*TEHPI7I«TEMP(7)*TEHM&I/(V«V*II*-ETxI) 
OFPOTaFM«TEMPI Ml«( l«-FHi*TEHp||2l 

TJOaTJ 

TJ«TJO-FF/OFFOT 

IFITJ  .LT*  3wO«l  TJ>300« 

IF  UBSI ITJ-TJOi/TJO)  «LE«  l.E-Hl  60  TO  76 
|TERTj>iTERTJ«i 
IF  IITERTJ  «6T.  201  6U  TO  733 
60  TO  73 
733  TJOaTj 
ITERTjaC 
TJa*S«(TUPaTLN) 

IF  (TJ  .LT*  300*1  TJB300* 

IF  UBSI ITj-TjOI/TJOl  *LE*  l.E-Bi  60  TO  76 
ITH>ITH«1 

IF  IITN  *6T*  201  CALL  SPRINT  UHE0S3  .AHERRUR3I 
60  TO  73 
C 

76  TL-TJ 
TT»TJ 

60  TO  f2 
C 

C •••  LIQUID  RE6I0N  ••• 

C 

7<l  1«FLA«IJI>32 
ITERTaO 

XXaTOLO/TEHPI 121 

71  F*TEMPI23UXX«ITEMPI10UXX*TEMPI22n«TF.MP{|3)*lAL06ULF*XXa|*0) 

« -ALF*XXI 

0F0XaTEMPII9U2*«TEHPI22l«XX-»TEMP||3l«ULF/IALF*XX*l*0l-ALF) 

XXOaXX 

XXaXXO*F/OFnX 

IF  IARSHXX-XXUI/XXOI  «lE*  t*E*RI  60  TO  72 
ITERTalTERTal 

IF  IITERT  *6T*  291  CALL  SPRINT  I4ME0S3  .AHLMRONII 
60  TO  71 
C 

72  TL«XX«TEMP| 121 

77  PCSa0*6«CC0NI2tNI«TL««2*(ESTC0NI |«*N>-«ANNAS I •TlNPI 2 1 

PRC  aPCS-TEMP 1 2 1 aTEPP 1121 aECONI 3 .N I • I TEMP  1 1 1 1 ♦6AMNA5 1 -TEMP  I * I 

• -TEMPl2l«TeMP(  12l«l  rCHP||<l)««AM»<aSl*IALU6IALFvXXa|*|-|ALF*l 

• AL06lALFa|«l-ALF«IXx-l*llaTrNP( I IFALF 
PIJUTEMPI  |«l*PRC«ECON(  /*N| 

TtUUTL 
60  TO  100 
C 

c •••  vapor  RTnION  ••• 
c 
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/ 


9C  IVFLA6(JI>3S 

TEMP(t»«-TEMP(2)*CC0N|l/»N) 

VS  lTERT-0 

TEHP( U)«-E( J)*TEHP(V) 

TEHP(25»«TEMP(2)*ECON<&,N) 

TEMP(  17)«EC0iJ(S,N)*TE(tP<25» 

TEHP(  18)«.5*TEMPU7I«TEHP<25> 

TT-TOLO 
ITH-C 
TUP>I .E6 
TLHO30S* 

VI  IP  (TT  .GT.  ECONIIV.NH  GO  TO  VII 

ve«ECON( I8|N J-rT*(ECO»M  lV,NI+EC0U(2CtN)*TT) 

TEMPI  7»«-EC0N(  lV,NI-2.*t;CON(20»»J)*TT 

TENP(8l>-2**kCON(2CtM 

60  TO  VI 2 

VII  V8«EC0NU,NI/(TT*«.2S) 

TEMP(7)«-.25*rC0NU,M)/|TT**U2S) 

TEHPI8l*S.«EajN(6.N)/l lo.«TT**2.2S) 

V|2  ETA-VB/V 

TENP(S)a(2.C-ETAI/< 1 1 .0-ETA ) aad I 

EV-TEMPII l•TT-2.•CC0N^S,N)•TT•TT•TEMPI5)•TErtPI7)/V  ♦ ESTCONIIOtN) 

. •! .E-li^IEMPIVI 

F-EV-EIJI 

IF  IF  .GT.  J.l  TUP«AHIN1 ITUP.TTI 
IF  IF  .LT.  0.1  TLW>AIIAX1ITL».TT) 

IF  I ITH  .GT.  SI  60  TO  V|3 

OF  -TEMP  I I |•H.•EC0NIS,^H•TT•7EMPI5)•TEMP»7)/V 

. -2.«eC0NI5,ni*TT*TT*TEMPI5)«IEMPIB)/^ 

. ♦2.*EC0NI5.HI*TT*TT*TEMPI7)»TEMPI7)/IV*V»I I .-ETa)«»3i 

. -A.*EC0NI5,NI*TT*TT*TEMPI7)»TE«PI7)«TEMpl5)/IV*V*l I.-ETAI  I 

TTO-TT 
TT«TTO-F/OF 

IF  ITT  .LT.  303.1  TT«30w» 

IF  I ABSIITT-TTOI/TTOI  .uE.  I.E-3)  CO  TO  V2 
|TERT«IT£RT+I 

IF  IITERTJ  .GT.  20)  60  TO  V.j 
60  TO  VI 
V|3  TTO-TT 
ITCRTJ-C 
TT".5«ITUP>TL*) 

IF  ITT  .LT.  330.)  TT"30C. 

IF  I ABSI ITT-TTOI/TTO)  .LE.  I.C-V)  60  TO  V2 
ITH«ITH*I 

IF  I ITH  .GT.  201  CALL  SPRINT  I6NE0S3  .AHERRUR2) 

60  TO  VI 
V2  TCMPI5)"ETA 

PV  •CC0NI5.N)«TT*TEMri2)*l I.C*TEPPI5)*(I .3*TEMPC5I*I I .-TEMPIS) ) ) ) 
. /II I .O-TEHPISI l••3l♦TEHP^^)•TEMPI V) 

VR  IF  IIJOIN  .NE.  01  GO  TO  VA 
PIJ)"PV 
Tl  J)«TT 
60  TO  ICO 
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C •••  MIXED  PHASE  REGION  ••• 

C 

195  lVFLAG(J)a36 
1TER«!TER*! 

IF  (ITER  <GT.  31  CALL  SPRINT  (AHES03  i^HERRURII 
VXH|aVHN( I (N) 

VHM|aVHN( I .Nl 
EXH|>EVHN( 1 ,NI 
EVHlaEVHN( I ,NI 

EEqMI>EVH|«(EVHX(3.NI-EVMN(2.Nli*l V-VHMn/(VHXI2iN)-VHN(2.NI ) 

C 

DO  197  JJ>2.20 
JM«JJ 

EEOaEVHN(  JJ.NIMEVHXC  JJ,N)>EVMN(  JJ.Nn«(  V-VHNIJJiNI  1/ 

. IVMX(  JJiN)-VHNUJ.NI  i 

IF  IE(JI  .LT.  EECil  60  To  I9(> 

EBH«tVMN(  JJ,NI  + IEVH1-EVHNIJJ,NI  ) • ( V- VMN  IJ  J , H I I / ( VMM  1 - VMi(  C J J tN  I I 
IF  (E(JI  .6E.  EFHI  60  To  32 

EBLaEVHXUJiNI^ILVMXCOJtNI-EXMl  ) • < V>  VHX  IJ  J . N I ) / I MX  ( J J I -V  XH I I 
IF  (E(JI-FBL)  199,32.32 

196  EEQMIaEEQ 
VXMl-VMXI JJ,M) 

VHMIaVMNi JJ.NI 
EXMlaeVMXI JJ.NI 
EVM|>EVMN( JJ.N) 

IF  (ABS(  VMXI  JJ-»|  ,NI>VIIN|  JJ«|  .Nil  »LE*  C»i  GO  TO  |9B 

197  CONTINUE 
C 

I9fl  EBLaEVMXI JM,NI*VMX( JM,N)/V 
IF  (EIJI  tOE*  EBLI  60  To  32 
C 

SLENa.EVMNI  JM,N)/(  I • /L'CoN  ( | 6 tN  | • 1 • / VMN  I JM  .N  ) ) 

E8HaEVMN(  JM,NI>SLENa(  I .C/V-1  »0/VMN  I Ji1  ,U ) ) 

IF  (E(J)  .6E*  EBHI  60  TO  3S 
C 

TEMPI  I laTMPI JM-2.NI 
TEMP(2>aTMP(JM-| .Ml 

TEMP(3)a(E( J)/EEOHI )a(TtMP(2l-3rC. »♦30C• 

GO  TO  101 
C 

199  TEMPI  I laTMPI JM-I .Nl 
TEMPI2laTMPI JM.NI 
TEMPINIalEEQMI-EI  jn/IEI  J)-EEOI 
TEMPI3lalTEMPI2l  aTEMPIRi^TEPPI  I I I/I  I .O^TEMPMI  ) 

C 

101  TEMPISIaALOSIPKESI JM-I.nI/PRESI JM.NI I/I  1.0/ TEMP  I I l-l ./TEMPI2I I 
PI JIaPKESI JM.NI aEXPITEMpiSiai l,o/rEMPI3l-I.C/TEItP(2l I I 
IF  (FIJI  .LE.  ESTC0NI2.HI  /i•E♦l2l  PlJiaO. 

TIJ|aTEMPI3l 
GO  TO  ICO 
C 

C PIJiajOIN  PRESSURE 

C 

96  XVa| .0«V«ESTC0N( I ,NI 
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TCHP( |S)-eSTCON( I .NI/( I.O-XV) 

TCMP(  |&)«XV*TeHP(|S)*(  I.E-6*ESrC0M(««NI/(  |.0>eSTC0N(  Uf.n*xvn»«2 
TEMP(17I«1 .S-XV^C I.3*.5»ESTC0N( I7iN»-.5*XV*ESTCON(  18tN) ) 
6AMMXS«ESTC0N( I7tH»-XV*tSTCON( I8|N) 

PI V«TEHP( I6)«TtMP( 17)4TfMPC I S I •CAHM AS* ( E (J I -ESTCON ( 22 . N I) 
PCSV*,5*EC0N(2iN)*TL**2*IESTC0»^(1  9 |N» -GAMMAS)  •TEMPOS) 

XXTM«TL/TMV 

TEMPI |8)*2.*( XV*ESTCOri( 16.N)-ESTC0N( 17 (Nl^.iJSS) 
PCTMV*PCSV-TMV*EC0N(3,)))*TEMP(  I S ) • ( TEMP  ( 1 8 ) *G  AMM  AS  ) 

• -.5*EC0N(H,M)*ECU)4VV-TEMP(1S)*THV*(TEMP(  18)*GAMMaS)*TEMP(  n* 

• ( At06(  Atr*XXTt1*l  .C)-ECON(  |S«t))-ALF*(XXTM>l  •:)  )/ALF 
PL*P1V*PCTMV*EC0N( 7,N) 

E0V|*EC0N(  1 1 .))) 

T0*30C* 

TTERMI-3. 

TTERM2«5» 

IF  (TT  .LE.  Cl. I SO  TO  7h3 
TTERMI«TT*ALOg(TT) 

TTERM2*TT*ALUGIABS(TT/To) ) 

783  CONTINUE 

TERMI«2.*TT*TEMP( I )*.S*tC0N(2,N)*TT*TT-2.*TEMP( 1 )*TTERM2 
. -ECON(2iN)*TT*(TT-TO) 

TERM2»2.*TT*TEMP( I )• (ESTCONI I7.N|-ESTC0N( IB,M) )+.S*ECON(2.N)*TT* 

• TT*ESTCON( I9.N) 

ASVI*COVI*TFKMI*2**TT*TlHP( I )*CSTCUN( I0iN)*( I .C-V I *ESTC0N ( Id))) 

. -TERM2*AL06< VI*ESTCON( 1 ,N) ) 

AS»E0V*TERMl«2«*IT*TEfiPt I )*ESTfOM(IBiN)*XV-TtRM2*ALOGi I .0-XV) 
AV*TT*ECOH(5iN)*  AL0G(  TtMM2  ) )♦!  7.-ETA*  ( A.*ETA  ) ) • . S/ ( ( I . -ET A ) • *2  ) ) 
. -TEMPI  I )*TTERMI-CCON«  I7(N)/V*ESTCON(  I3.N)*I  *c:-12 
ETAVI*VB/VI 

AVVI*TT*CCONISiM)*I ALOGt I */V 1 ) * I 7 .-E T A V 1 • I A . *ETm V I ) ) • *B/ 

. I I I .-ETAVI )**2) )-TEMm( I I^TTEFNI-ECONI I7,N)/V) 

• ♦ESTCONI Ii,N)*I.E-l2 
XXTM*TL/EC0Ni8tN) 

AUVI*ASVl4Cr0NI IO.N)-«S*ECON(9«N)*ECONI3iN)-ITL-ECONI8iN)  )• 

• IEC0N(3.N)^TEMP( I ) )*2»*TEMPI I )*ECONr8.N)*I ( ALF*XXT(KI  .)• 

• ALOGI ALF*XXTM4I. |-IALF^I,I*At06I ALF^I. ) )/AtF 
XX*TL/TMV 

ECI0V*EC0NI8iNl*»5*CC0NVV 

AL*AS^ECI0V-»5*EC0N9V  *ECON( 3«N )- I TE-TMV ) • 

. IECON(3iN)^TEMPI I ) )♦^.•TEMO| I )• TMV*I  IALF*XX^I  . )• 

• ALOGI ALF  *XX^I • )-l ALF^I .)*AtOG( ALF^I . ) I7ALF 
P(  J)«FM*PV^I  1 «C-FM)*PL^I>FMOV*(AI-AV-ALVOAV91  ) 


TI J)*TL 

ESC222? 

100  EIJI*I*E^I2*EU)  ESC2230 

PI J)*I »E^I2*PI J)  ESC228D 

RETURN  ESC22S0 

END  ESC22AC 


r%  r% 


APPENDIX  C 


LISTING  OF  DATA  GENERATION  ROUTINES 


SUBROUTINE  OATSCN 


COMMON 

AMASS! ICI .OELTMI (Id .OELTMNI 131 •HaTBNUIICI tHATLlO( IO> • 

lEO 

30 

• 

RATIO!  Id  tNMTRI.S 

lED 

SC 

COMMON 

VlSCOL!  101  f VISCOvI!  131  , COItRST  , OELT  1 M , M A XC  YL  , MOOE , NRESRT  t TMA  X 

lEO 

S3 

COMMON 

(PLATE, UZERO 

lEO 

<•3 

COMMON 

ALPHA !8I 1 ,AMU( I 1 I ,RBT£HP!2) ,CUNVRC • ENERSY ( 2 1 • MAXPAS , NALPH A • 

lEO 

7C 

• 

NFLUO!  lai  ,NPE!  ICI  .NSCTTUd  ,SSTART!2)  ,SST0P(2)  .NONITI  , 

lEO 

ac 

• 

NUN|T2,NUN|T3,NUN|TS,NriT| li.  I , NL06E . H 1 NC , NKF 1 LT , NFUNCT 

lEO 

9? 

C IMUN 

STIMEI ,ST IME2,CNRGY,niOTH 

lEO 

I3C 

COMMON 

I0UNIT!6) 

lEO 

1 1C 

COMMON 

NCYEPT  !2CI  ,TE0|T!2r 1 .NABSISI .rtOROIS) , i ZPLOT ( S I ,LPLOT  ! 5 I • 

lEO 

123 

• 

IPOT  ,H0MEt4T,HDUMH,iiFE0IT,NPRlNT  ,TPR  I NT  ,NsEO  1 T 

lEO 

I3C 

COMMON 

ICUTL,  ICUrR,lAOOl.,lAOOR,lDDEM,Oi!ZI  .OKATIO.AOOMaS.OUMI 

lEO 

ISO 

COMMON 

I VZOMEtVWIOTH 

lEO 

ISC 

COMMON 

OELTAM! ic:0l .OIScPT! 121 .MERRORIZ) , ITEMPI25I .MZONE! IO,SI , 

lEO 

I6C 

• 

TEMPI2SI ,0TAVG,0TN,0TNPI .ICOURN. JSAVE,MAPFL6,MAXZ,MAXZP1 , 

lEO 

|7C 

• 

NCYCLt,NSAVE,POTpOS,POTM».G,TIME 

lEO 

IOC 

common 

C!  IC.^dl  ,E!  10C3I  .hRATE!  ICt3,2)  tlBFLAG!  lOQC  I ,1  VFLAGI  lOOOl  , 

lEO 

190 

• 

OLOHHO! ICCCI ,P! I„33l .0 ( I'C? I ,RHO I ICOd ,S ( 1 OOC I ,SU( I 003) , 

lEO 

2C? 

• 

U(  lOSJI  fXdOCSI  , BLANK  !39| 

lEO 

210 

lEO 

22C 

COMMON 

/EOST/  ESTC0N(33,3CI 

COMMON 

/TEM/  TUCCJI 

COMMON 

/TSAVE/  TOLnJtCOtOJ 

OIMENS 

ION  rHP(23l  ,VMX(23)  , VKN(2*:  1 tCVMX  ( 20 1 ,E  VMM  ! 20 1 , PRES!  20  > 

• ••  calculation  or  critical  poi».t  ••• 

c 

JAA«I I 

JABaJAA^I 

AJAAaJAA 

NKat 

NK|aNK-l 

NK2aNK>2 

NK3aNK*3 

NaNSAVC 

jSAVEal 

JSaJSAVE 

VOSal ./CSTCONI I iN) 

TH0aESTC0N(2l ,NI 

0ELTa2,»TM0 

TJSaOELT 

TOLOJaTjS 

0ELVa3.aV0S/20* 

ViNaVOS 

VCESTa<t.S«VOS 

IR-I 

ITa| 

112  ITalT^I 

RHOI JSI>I*/VIN 
VJSaVIN 
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^hkbd 


CALL  ET  ( VJS.TJStCJS) 

E( JSi-EjS 
CALL  E0S3 
PI »P« JS) 

JVC»I VCEST-VJSI/UELV 
00  110  J«l t JVC 
VJS»VJS40ELV 
RHOI  JSiaU/VJS 
CALL  ET  (VJS.TJS.EJSI 
E( JSi-EJS 
CALL  E0S3 
P2»P« JSI 

IFHP2«Pn/0ELV  .GE.  C.l  GO  TO  III 
I 13  PI-P2 
IK-IK^I 

IF(IR‘>2  *GT.  C I GO  TO  l23 
DELV«0.7*OELV 
GO  TO  112 
123  DELTa.5«0ELT 
TJS-TJS-OELT 
TOLOJaTjS 
IR«I 

GO  TO  112 

III  IFIIT-ICO  .GE.  0 I GO  Tu  IIS 
TJS»TJS*.5«0ELT 
TOLOJaTJS 
VIN>VJS-2««PELV 
GO  TO  112 
IIS  TCaTJS 
VC>VJS 

CALL  ET  < VJS.TJS.CJSI 

ECCaEJS 

PCaP2 

PRINT  2C0iTCiVC.PCiECC 
C 

C •••  COEXISTENCE  CALCULATION 

C 

00  100  ja2iJAB 
AJaj-l 

TJSaTC-AJalTC-THOI /AJAA 
TOLOJaTJS 

OELRHOa.CI aCSTCONI I ,N| 

RHO( JSIaO* 

PlaO. 

K|a| 

00  20  Kaliino 

RHO( JSIaRHOt JSI«OCLRHO 

VJSal ./RHOI JSI 

CALL  ET  ( VJSiTjSiEJSI 

E( JSIaEjS 

CALL  E0S3 

P2aP| JSI 

IFIIP2  (LE*  Pll  (AND*  Ul  .EO*  211  GO  TO  20 
IFIP2  *GT.  Pll  GO  TO  |2 
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PHAXaPi 
Kl>2 

GO  TO  2C 

IFIKI  (EQ*  l>  GO  TO  20 
PNIN^P2 
60  TO  U 
2C  PI«P2 

16  PE0«*S*IPHAX4>PNIN) 

IMPEQ  *6T.  &•!  GO  TO  22 
PEQ-.S*PHAX 
22  00  25  KIT>I ,50 
00  330  K-l ,NK2 
00  *«CC  L»l.« 

AL>L 

VJSX |C***(NK-<) l/AL 

RH0(J)«l «/VJS 

CALL  ET  IVJSiTJS.CJSI 

EIJSi>EjS 

CALL  E0S3 

PI>PI JSi 

IFIPt  •GE«  AEQ)  go  to  Rl 
SCO  CONTINUE 
300  CONTINUE 
K>NKI 

00  RS  L«l3i2S30 
AL>L 

VJS>IC3./AL 
RH0IJSI>I*/VJS 
CALL  ET  IVJStTJS.CJSI 
El JSI>EJS 

CALL  E0S3 
PI«PIJS) 

IFIPI  .GE*  PEQI  GO  TO  Rl 
RS  CONTINUE 
Rl  RH0M|N«RN0( JSi 
SUN2>0* 

IN*  I 

IFIK  «GT*  NK3I  GO  TO  S3 
00  55  KK*K,NK3 
OELRhO*I«/I  IO.««iriK-KKi  ) 

00  6C  LL*LflO 
AL*LL 

RHOIJSi*AL«OELRHO 
VJS*l«/RHOI JSI 

CALL  ET  I VJ5,TJS,EJSI 
El  JSI*EJS 

CALL  E0S3 
PI>PI JSI 
F«P|-PEO 

SUN2>SUH2«F«UELRHO/RHOI JSI*«2 
IFIF  .LE*  9*i  IN«2 
60  CONTINUE 
55  L>l 

S3  DO  65  LL>L,2530 
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OELRHOa.Ol 

AL-LL 

RHO( JS)>AL/ICa* 

VJSal ./RHO( JS) 

CALL  ET  ( VJS.TJS.EJSJ 

E( JSI-EJS 
CALL  EOS3 
PlaP( JS) 

Papi-pEO 

SUR2aSUR2>F*0ELRH0/RH0( jS»a»2 


1F(  IN 

• EQ» 

2) 

GO 

TO 

69 

1F«F 

• GE* 

0.) 

GO 

TO 

65 

lNa2 

6H 

1F(F 

.GE* 

OtI 

GO 

TO 

62 

65 

CONTI 

NUE 

62 

RHONAXaRHO 

( js 

I 

1F(K1T  .GEt  2)  GO  TO  A8 
PEQlaPEQ 

1F(SUH2  .LTt  C*l  GO  TO  72 
PE02aPEQl^< l*(PMAX-PEOI ) 

71  PEQaPEQ2 
SUMlaSUR2 
GO  TO  2S 

72  1F(PN1N  .LE.  0*1  PHlNaa. 

PE«2aPEQI-. I«(PE01-PM1N) 

GO  TO  71 

«8  IF  ( ABS(SUHt-SUM2l-l*C>6  tLE*  0*1  GO  TO  69 
PEQ3a(SUM2*oC01/SUMl-PEa2l/ISUM2/SUhl-l . ) 
1F(PEQ3  *LT.  C*l  PCQ3a.>iapEG2 
PEQaPEQ3 
PEQlaPE02 
PEQ2aPEU3 
SUHlaSUN2 
25  CONTINUE 
69  VNN( J|a|./RHUNAX 
VhaVNNl Jl 
VNX( J|al  ./RHONIN 
VXaVNXUI 
PRES( JlaPEQ 
TNP( JlaTJS 
CALL  ETI VN.TJS.EJSI 
EVNN( JIaEJS 
ENaEVNNl J I 
call  ET(VXiTJStCJSI 
EVNX( JlaEJS 
EXaEVNXl Jl 

IF(PE0-1«E-S  (LE*  0*1  GO  TO  9 
IOC  PRINT  5CCtTJ:i,VN.VX,CMiEXiPC0 
JaJAB 
9 JMAXaJ 

PRINT  SCOtTJStVNiVX.CMiEXiPEN 
VNX( I lave 
VMN( I )aVC 
EVNX ( I laECC 
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EVHN(II>ECC 
PRES( 1 l»PC 
TMPC  I »-TC 

200  F0RMAT(5X,5HTC  > •EIO.S.SX.SHVC  • tElO«5.5X  .bHPC  > lEIO.Si 
25X.SHEC  > . .ns. 5. ///I 

500  F0RMAT(5X.6HTHP  a . FI  3. 5 »5X .AH VMN  a ,£ 1 0 .5 . 5X  . AH VMX  > .EI0.5.SX, 
27HEVMN  a .EI0.5.SX.7HCVMX  a .CI'.S.SX  .7HPRES  a ,E19.5./) 

RETURN 

END 


i 

I 


c 


SUBROUTINE  ENCALC 


I 


C 

c 

c 

c 

c 


I 

! 

c 

c 

c 


ENC  10 
ENC  2C 


COHHON 

COHHON 


ARASS(ICItOELTItlllS).OELTNNI|jltMATBNU(i:iliriATLIOllOli  ENC  3C 

RATIOI  lOi  (NKTHLS  ENC  ‘<0 

ViSCOLt  101  .ViSCOUl  tai  .C0llRST.0ELT|*4,MMXCYL.N00EiNKESRTtTNAX  ENC  SC 


COHHUN  IPLATEiUZERO 

CONhON  ALPHAtni I iANUI 1 1 1 iBBTENPIZt  fCUNVRG.EuERGY(2l tNAXPASiNALPHA. 
. NFLUOl iul iNPEI IGi •NSCTTIICI .SSTARTI2) .SST0P(2I iNUNITl . 

. NUNIT2tNUNIT3itiUN|TH,NFITI  Ul  .NEOGr.NINCiNKFILT.NFUNCT 

COMMON  STIMEI  •ST|l1E2iCNi<GYtli»IDTH 
COMMON  IOUNITIAI 

COMMON  NCYEnT(20l .TE0IT(23t .NARSISI tNOROiSi tlZPLOTISI  iLPLOTtSI t 

* lPOT,NOMENT,NOUMH,NFEOIT*NPRIMT,TPRINT,HaEOIT 
COMMON  ICUTL.ICUTR.IAOOLiUrOR.IOOEN.OMZI .OH XT  I 0 . ADDMAS ,OUMI 
COMMON  I VZONE.VIKIOTH 

COMMON  DELTAMI ICjOI .OISCPTI 121 .MCRKOR I 2 I . I TcHP ( 25 i .MZONE (10 iH I . 

. TEMP  1 25  I .OTAVG.OTN.OTNPI . ICOURN . JSA VE .M APFLG .HaXZ ^MAXZP I • 

* NCYCLE.NSAVE.POTHOS.POTNfG.TIME 

COMMON  C(  ia:OI  IEMC30I  .HRATE(i:';.2t  .IGFLAGdOCti.lVFLAGlIOOOl  . 

* OLORHOI lOCCi .P(iu3St .olircoi .HHOI lOOOl .S( lOOOl .S0( 1000) • 

. U(l0CtJliX(l500).ULANK(39) 

COMMON  /CEOIT/  OAMAGEI Io30) .NSPALLtlO) .SKAPPAl 1C) .SLAMOA( 10) . 

* SINITLI 10) .SFINALi 10) .ICOMHZ.)SPALL 

COMMON  /EQST/  ESrCON(30.30) 

COMMON  /TEM/  T(ICCO) 


ENC  *? 
ENC  73 
ENC  BO 
ENC  VO 
ENC  103 
ENC  DO 
ENC  120 
ENC  130 
ENC  )‘<C 
ENC  ISO 
ENC  |t>0 
ENC  170 
ENC  tdO 
ENC  IVO 
ENC  233 
ENC  210 
ENC  220 
ENC  230 
ENC  2B0 
ENC  250 
ENC  260 
ENC  27C 


COMMON  /TSAVE/  TOLOJ.COLOJ 


SOTl-C* 

S0T2«0* 

SMASSaO* 

FRAC«C> 

FRAClal* 

FRAC2>1. 

•••  CALCULATE  SOURCE  TIME  STEP  ••• 

IF  (NOOE  «EQ«  I)  GO  TO  20 
TEMPI) )»T|ME-0TNP| 

IF  (TEMPO)  «6E.  SSTOPOD  GO  TO  13 
IF  (TIME  «LE*  SSTARTOIi  60  TO  1C 
SOT )>0TNP) 

IF  (TIME  .GT.  SSTOPO)  .AND.  TEMPO) 
. TEMPO) 

IF  (TIME  «GT.  SSTART(I)  .AND.  TEMPI  I 


ENC  2B3 
ENC  290 
ENC  300 
ENC  >10 
ENC  32C 
ENC  33C 
ENC  390 
EMC  3SC 
ENC  3*0 
ENC  373 
ENC  3BC 
ENC  390 
ENC  900 
ENC  910 

.LT.  SSTOPOD  SOT|>SSTOP(  I )•  ENC  920 

ENC  930 

.LT.  SSTAkTOII  S0TI>T|ME-  EMC  990 


• SSTARTO) 

1C  IF  (TEMPO)  .GE.  SSTOPIZD  60  TO  23 
IF  (TIME  .LE.  SSTART(2I|  60  TO  20 
S0T2»0TrjP| 


ENC  950 
ENC  9*0 
ENC  970 
ENC  9B0 


IF  (TIME  .GT.  SST0PI2)  .AND.  TEMPO)  .LT.  SST0P(2)I  S0T2«SST0P(2)>  ENC  99; 

TEMPO)  ENC  SOO 
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I 


rt  rt  n f\  rs  rs  r%  rt  r>  n n 


IF  CTiHt  .6T.  SSTARTI2»  .AMO.  TEMPIII  .tT.  SSTARTIZII  SoT2«T|ME- 
. SSTARTI2I 

•••  CALCULATE  internal  ENER6T  ••• 

20  M>| 

N>MATLiDIHi 

IOONE>C 

lEOSaESTCONI3C.NI 

IF  IIhOOE  .NE.  0 .OR.  (mode  .EO.  0 .AND.  NMTKLS  .EO.  Ill 
. .ANO.  TiME-OTNPl  .6E.  ANAXI ISSTOPI 1 ) .SST0PI2) I)  lOONEal 

•••  Main  oo  loop  for  calculating  external  heat  sources,  internal 

ENERGY.  UEWUTORIC  STRESS  COMPONENT.  ARTIFICIAL  VISCOSITY 
COMPONENT.  PRESSURE.  AND  STRESS.  ••• 

IF  (NFUNCT  .nE.  21  GO  Tu  30 
FRACIabOUNUmHE.il 
FRAC2-BOUNOITIME.2I 
30  00  160  ja| .MAXZ 

IF  tj  .lt.  matbnoimii  go  to  no 

HaM«i 

NaMATLIOlHI 

IEOSaESTCONI30.NI 

SO  IF  (ABSIUIJII  .GT.  0.1  GO  TO  60 
IF  IIOONE  .EG.  01  GO  TO  60 
Jlaj 

00  SO  KaJI ,HAXZ 

IF  lApSlUIKII  .GT.  e.l  GO  TO  60 
SC  CONTINUE 


ENC  SIO 
ENC  S20 
ENC  S30 
ENC  SSS 
ENC  SS7 
ENC  S60 
ENC  S70 
ENC  seo 
ENC  SVC 
ENC  600 
ENC  6IC 
ENC  620 
ENC  630 
ENC  6S0 
ENC  6SC 
ENC  660 
ENC  670 
ENC  680 
ENC  690 
ENC  700 
ENC  710 
ENC  720 

ENC  730 
EnC  7SC 
ENC  7S0 
ENC  760 
ENC  770 
ENC  780 
ENC  793 
ENC  8CC 


GO  TO  I7C 

60  IF  IIBFLAGIJI  .E0»  I .OR.  IBFLAGIJI  .EU.  31  GO  TO  160 
NSAVEaN 
JSAVEaJ 
EA0D«0« 

IF  IHOOEi  70.60. 90 

70  IF  IIOONE  .EG.  I .AND.  NCYCLE  .NE.  II  60  TO  90 
SHASS«SMASS6DELTAMIU| 

CALCULATE  HEAT  ADDITION  FROM  ELECTRON  BEAM  SOURCE 

CALL  TBEAH  (SHASS.SOTI .FRaCI 
EAOOaFRACtHRATEU.I  I 
60  TO  90 
C 

c CALCULATE  HEAT  AOOlTlON  FROM  RADIATION  DEPOSITION  SOURCES 

C 

80  EAOOaSOTlaHRATEtJ.I |arRAC16S0T2»HRATEIJ.2l#FMAC2 
90  OUaUtU«| l-Ut Jl 
ITERaO 

OVaOUaOTNPI  /OELTAIttJI 
EIJIaANAXI lEI JI«EAOO*.SaSI JI6DV.0.I 
C 

C CALCULATE  0EV|ATUR|C  STRESS  COMPONENT  ANO  ARTIFICIAL 


ENC  8IC 
ENC  829 
ENC  83C 
ENC  819 
ENC  8bC 
ENC  e»s 
ENC  870 
ENC  880 
ENC  890 
ENC  900 
ENC  910 
ENC  920 
ENC  933 
ENC  910 
ENC  9S9 
ENC  960 
ENC  97C 
ENC  980 
ENC  990 
ENC 1 009 
ENCICIC 
ENC 1 020 
ENCI030 
ENC 1 010 
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n n r>  o o o no 


VISCOSITY  PRESSURE  COMPONENT 

IP  (ESTCON(StN)  .CT.  C«)  CALL  STRESS 
IF  (ESTCONietN)  .LT.  0*1  CALL  STRCSl 
Of JIaQFUNCTIN) 

CALCULATE  PRESSURE.  STRESS*  AND  INTERNAL  ENERGY 
Ef JI>AHAXIfCf JI*(.S«SD( J)«QI Jt t*DV.C*l 

•••  FOR  OATGEN  RUNS.  REMOVE  C FROM  THE  NEXT  THREE  CARDS  ••• 

El J)>0* 

CALL  OATGEN 
CALL  EXIT 

GO  TO  mO*ICO.IC3.IIO).  lEOS 
IOC  HOLOP-PfJ) 

El JI>AMaXI (CIJ)-*S«PIJ)*OV.O.t 
110  CALL  EOS 

Sf  J)>PU)«SDIJ) 

GO  TO  (I23.I30*I3C.I20).  lEOS 
120  EU)>AMAXI  (El  J)-*5*P(  J)*0V*3*I 
GO  TO  ISO 

130  EI>AMAXI (E(JI-«S*(P(J)*HOLOPt*OV*Q*) 

IF  (ABSI (Cl-ei J) I/El ) .LT.  S.E*3I  GO  TO  IVQ 

ITER>ITER«I 


ERROR  EXIT 

ERROR  I « ENERGY  ITERATION  HAS  FAILED  TO  CONVERGE 
IF  (ITER  *6E*  10)  CALL  SPRINT  (AHENCALCiAhenHORM . . 


C(Ji>CI 
HOLOPaPI J) 

60  TO  no 

MO  CONTINUE 

SET  momentum  editing  FLAGS 


TEMPI  I l•TINE-0TNP| 

IF  ITeMPII)  «6E*  SSTOPID  .AND*  TEMPII)  *GE*  SST0P(2)I  GO  TO  ISO 
IF  (EIJ)  *GE*  ESTCONIIJ.NI)  MIONUM.IMJ 

IF  (EIJ)  .GC*  ESTCONIII.N))  HZ0NC(H,2)>J 

IF  (EIJ)  •GE*  ESTC0N(|2,N))  MZ0NEIM,3)bJ 

IF  (EIJ)  tGE*  ESTC0NII3.N))  NZONCIM.Daj 

CHECK  FOR  SPALLATION 


I ■ 4- 


CNC  IS2C 

I5C  IF  (ISPALL  •NE*  0)  CALL  SPALL  ENCIS3C 

l&C  CONTINUE  ENCISH? 

I7C  return  ENCISSO 

ENCIS6C 
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r%  f%  r\ 


c 


SU8R0UTINC  Cr  I VJS«TJS«L'JSt 


C 

C 


C 

c 

c 

c 

c 


ESC 

COMMON  AMASSI tCI tOELTMl I I C t lOELTHN 1 1 3 1 .MATBNU I I S I iMATL I D I 101 . CSC 

. MATtOI tot iMHTRLS  CSC 

COMMON  VtSCOLlIC) fVISCOwl 101 iCOURST, belt IMiMaXCYL (MODE •(iMCSMT.TM AX  CSC 

COMMON  tPLATC.UZERO  CSC 

COMMON  ALPHAien  tAKUl 1 1 i . BBTEMP 1 2 1 tCUNVRG , CNLRG V 1 2 1 «M A XPAS .NaLPHA . CSC 

. NrLUOIIO).llPC(IOJtNSCTTII31.SSTART|2)<SSTOPl2).NUNlTI.  CSC 

• NUNtT2.MUNIT3tliUNlTN,NFIT(ICIiNEDGEiNtNCiNKFILT.:4FUNCT  CSC 

COMMON  STIMEI .STIME2iCNkGY,«|0TM  CSC 

COMMON  lOUNITU)  CSC 

COMMON  NCYCnTI2CI .TED I T ( 2C I »NARS ISI . NORO (S 1 . t ZPlOT I S I .LPLOT I S I • CSC 

. tPOT, moment, NOUMe,NFEOlT,NPRlNT.TPR|NT,NSEOlT  CSC 

COMMON  ICUTLtlCUTR.lAOOL, lAOOR.IOOCM.OMZl .OKATIO.AnONAS.OUMI  CSC 

COMMON  |VZ0NL,V*I0TH  CSC 

COMMON  OCLT AMI  I C 3C 1 , 0 1 ScPT 1 1 2 I ,NERROH I 2 1 , I TEMP ( 2S 1 .MZONCI IO,NI  • CSC 

. TCMPI2S),OTAVG,OTN,DTNPI,ICOUKN,JSAVk.,MAPFLG«MAXZ,MAXZPt,  CSC 

. NCYCtE,NSAVC,POTMOS,POTMtG.T|MC  CSC 

COMMON  Cl  IC*iOI  lEllCCCI  ,rlRATEI  i:*‘9.21  .IBFLAGI 10031 ,1  VFLAGI  lOCOI  • CSC 

• OLORHUIIOaOliPllUdCl*Qtir:COI,RHOIIC03l,SIIOOOI«SollOSdl«  CSC 

• Ul 10331 iXI lOOCl •nLANKt3«|  CSC 

ESC 

COMMON  /COST/  ESTCONiaC.aOt  CSC 

CSC 

COMMON  /GRAY/  T MP I 23,30 1 • VMX 1 2S*3C 1 • VMNI 20* JO  1 ,E VMI I 20 ,301 , 

. CVMNI23«301 ,PRESI2J,33I 


23 

3C 

NO 

SO 

bO 

83 

92 

ICO 

Its 

I2S 

130 

INC 

ISC 

163 

172 

IBS 

190 

20C 

210 

220 

230 

2NC 


DATA  GANMAF,ALF/|.0&S6i.  1/ 


OAT  A B01.TZ,AVN0,PI/I*3BL-I6,6.32C23,3*INI6/ 


COMMON  /CC/  CCQNI20i3CI 

3>3SAVC 
NaNSAVC 
RHOUIal  «/VJS 
r0L0>T3S 

IF  ICCONII.NI  .GT.  0*1  GO  TO  33 

•••  INITIALIZATION  OF  CWUATfON  OF  STaTC  CONSTANTS  ••• 

CCONII iNlalt/CSTCONI I ,N| 

CC0NI2iNI>CSTCUNI23,NI/LSTC0NI2N,NI 

CC0NI3iN)>CSTC0NI29,N|/tSTC0NI2N,NI 

CC0NIN,NI>I«I6JC’*S/CSTC0NI2N,NI 

LC0NI«,NI>0.|3NC-S/CSTCUNI2N,NI 

EC0NI7.Nl>N.SCN«CC0NI2iNmSTC0N||  ,NI*ICSTC0NI  l7fNI-CSTC0NI  |9,NI  1 
IF  ICSTCONI23,N)  .GC*  C«l  CSTCOMI 23 ,N I ■>9S0>* I B t I 3NL-S«S0«* 
CjTCONIZOiNI I/CSTC0NI2N.M 
CCONI  |2.,NI>t.-eSTC0NI|,Nl*C$TC0N(IN,NI 
EC0NI|S,NI>IALF«I(0I«ALUGIALF«I.CI 
TEMPI  I laCSTCONI |7,NI«|. J33 
TCMPI21>CSTC0NI2I  ,NI/I  l.3’*CC0NII2,NI  }**2 


ESC  Z63 
CSC  373 
CSC  380 


CSC  390 
CSC  NOC 
CSC  NIO 
ESC  N2C 
CSC  N33 
CSC  NNC 

CSC  NSO 

CSC  N7C 
ESC  NVC 
ESC  5CS 
CSC  513 
CXVI 
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TEn?l3lB|*«EC0NI |2,r«l«l2«*TCHP| t l♦ECON( |2tNI*ITeMP( 1 )• 

• I2*«ESTC0N( |7,NI-|.e667l«ESTC0M( 1 8 .N lit 

EC0NI«.NI>TCHP(2I*TEHP(3I  THI/I 

TEnPMIa2«*(ECCN(  l2.NI*ESTC0NCI**tNI*ESTC0N(  l7tNM«3333l 
TEnPISIa|.2*(  {•-ECONC  |2.NM*>*E-I2«ESTC0N(8.NI«*2 

EC0N(»,N»»EC0N(3.N»*(TEMP««H*TE»'P<2>*TEMP|3n**2/TEMP«5»  TI5VI 

TC*ESTC0N(2«»NI 

VBaCSTC0N(26tNI 

SI6HA>l3«*ESTC0N(2‘ttN|«VB/IP|*AVN0«8«) l••.3Ji333^3 
OT>(6**ESTCONl2<ltNI«VH/(P|«AVNOl  l•••3333333J 
EPSBB0LTZ«TC*(0T/l6AnnAF«SIGHAI 

EC0N(6tNlBPI*AVN0*ML.>:MHAF«SIGHA*(EPS/B0LTZl«**i:8333333l**3l/ 

• («*?*ESTC0N(2*I(NI  I 

FCON(  |C|N|s«M3*ECONIO.HI«*S«ECON(y  tU) 

TEMPI IIbESTCONI IA.NI«rCON( I2.NI 

TenPI2ls«5*l  l•E-A•ESTcOM(A,ttt•ECONn2•tnl••2/l  I.-TEMPM  I I 
TEnp(3l>U*TEMP(  I !•(  I.*.S«TEMPI  11*11  *-LSTCOi4 1 I 7 .N  I /ESTCON  1 1 6 .N  1 1 I 

• 73* 

ECONI 1 1 •Nl*TEHP(2l*TEtlP(3l*ESTCONI23*NI*l |.*ESTcONI 17|NI*  Tl IVI 

• ECONI  l2tNn*ESTCOtri22iNI 
TEMPI  251 -TEMPI  21 *CCON| 8, Nl 

ECONI  |7,NI-F5TC0NI27,tll/ESTC0MI2<I.NI**2  CSC  750 

IF  lABSlTMPI I «NI I .LE.  0*l  GO  TO  22 
DO  20  11-3.23 
JNAX-I l-l 

IF  iabsithpiii.nii  .lc.  o.i  go  to  21 

20  CONTINUE 

21  ECONI I A .N I -VHN UMAX. N I -PRES  I UMAX .N I • I VHN IJH A A- I . N I -VMN I JH AX .Nil/  I VMNOI 


• IPRESI JHAX-I .NI-PRESI JMAX.NI ) 

22  ECONI  |*i,N|a2.*ESTC0NI2l  ,NI 

ECONI  I8.NI-I.  |*EC0NI6.N|/IEC0NIM.NI**.25I  IVBOl 

ECONI |y.Ni-2.*EC0NI I8.NI/EC0NI I«I.NI-2.25*EC0NI4,NI/ 

• lECONI |R.Nl**l.2bl  tC2l 

EC0NI20.NI-|.2S*EC0NI«.nI/ICCON| |H.NI**2.25>>CCUNI  IB.NI/ 

. lECONI IR«NI**2I  IC3I 

PRINT  25.  IECONII.NI. 1-1. 201 
2B  FORMAT  I//.7H  ICONS .// .8 1 RX . |P5C IR*«./ I .// ) 

C CSC  «3S 

30  RNU-ESTCONI I .NI/RHOI Jl 

EX-I.*RnU  CSC  980 

C CSC  97C 

TEMPI  I laESTCONI |8.HI*CX  ESC  .993 

TCNPI2la.S*ll,  -8*CSTC0HI8.NI*EXI**2/| I .-TEMPI  I I I ESCIOOC 

TEMPI 3la| .*TEHPI I I • I I . • .S*TEMPI I I *1 1 .-ESTCON 1 1 7 ,N I /ESTCON I 18 .N I 1 1 / ESC  1010 
. 3.  ESCI023 

TEMPI  II .aTENPI2l*TEMP|3l*ESTC0NI23.NI*l l.*ESTCONI |7.NI*EXI*  ESC  1039 

• ESTC0NI22.<'M  ESCI?*:? 

IF  lEX  .GT.  0.1  GO  TO  8C  ESCI059 

TEMPI  I laESTCONI  I7.M-I.J333  ESCIC83 

TEMPl2laESTCUNI2l.NI/l l.-EXI*«2  ESC  1070 

TENPI3la|.*EX*l2.*TCHPI | l*EX«ITENPI I I«|2.*ESTC0NI I7.NI-I .88871-  ESC  1 089 

• ESTCONI I8.NI 1 1 ESCIC90 

GO  TO  53  ESCIlOO 

90  TEMPI  I |a2.*ICSTC0NI I7.NI-.3333I  ESCIIIO 


i 

f 


71 


non  n non 


TEMP(2(«ESTC0N(2I |N( 

TEMP(3(«|«-»EX*(TEMP(  l)«EX«nESTC0Nll7.N(-.33i3l«(TEMP(  I l«l*  (• 

• ESTCONI ISiNII I 

5C  TEMPJ I>J«TEMP(2I»TEMPJ3J 

TEMP(  m«2.«(EX«ESTC0fM  18,NI-ESTC0NI  I7.NI*.3333| 

TEMPI  I l•l•2•(  I .>EXI«(  .E>(2«ESTC0N<4,.>M««2 
TEMP(2(«i. 

IF  (EX  .GT.  u*(  TEMP(2(a|  ••»EX«(<f.«ESTCON(  )4»N(-(  • I 
TEMPI IS(aEC0N(3iN(«(TEMp( M|«TEHPri2) I **2/ I TEMP  I 1 1 *TEMP 1 2 1 1 
TEMPI  I laTEMPI (2I>.S«TCMP( IS! 

TMI«TEMP( I I 
V«(  •/RHOIJI 

IF  (V  *GT.  VII  60  TO  70 

IF  (TOLD  «GT.  TMK  GO  Tu  43 

•••  SOLID  REGION  ••• 

I VFLXGIJI«3Q 
TEMPI  I l■3.•EC0N(S.N) 

EJSaTEMPI II («T0LD«(*S«Ec0N(2tN(«T0LD«TEMP(  I ) ) 

GO  TO  (CO 

40  TEMPI  I (aTEMPI  I2I-».S«TEMP(  IS! 

TM2-TEMPI I I 

IF  (TOLD  tGT*  TM2)  GO  Tu  70 
•••  MELT  REGION  ••• 

(VFLAG(JI«3I 

TEMPI  I lalTOLO-TMl  I/ITI12-TMI  I 
TEMP(2(*EC0N(3.NI«EC0N(H«NI 
TEMP|3I»3.*EC0N|5,NI*TEmP| I (•TEMPI2I 
TEMPI5I«TEM*I I l••2 

EJS««5»ECONI2,N|*TOLD**2>TEMPI3I*TOLP*TEMPI ) |) -TEMP  1 5 I *TEMP  1 2 1 • 
. .S«TFMP(|SI 

GO  TO  ICS 

70  TEMPI  I laltSaECONISiNl 

TEMPI2I«ESTC0NI  I iNI/(l.k,<*EXi 
TEMPI  4 I ••  l‘♦3•EC0NI5tN)•,5•TEMPl  151 
TEMPI iai«2.  i^TEMr ; ! I 
TEMPI  131-TEMPI  I l/XLF 
TEMP(22|««S«EC0N(2.NI«TEMPI (21 
VI«ESTCON(  m.NI 
V2aESTC0N( ISiNI 
V>l .O/RHOI J( 

IF  IV  *LE«  VII  GO  TO  71 
IF  (V  .GE«  V2I  GO  TO  90 
C 

C •••  JOIN  REGION  ••• 

C 

IVFLAGI Jl>33 

TEMP(3I«2.L'«ECCN(S|N(«TEMP(2( 


ESC!  UO 
ESC! 133 
ESCI 110 
ESCI I5C 
ESC! I4C 
ESCII70 
ESCI lao 
ESCII9C 
ESCI23C 
ESCI2I0 


ESCI230 

ESCI320 

ESCI330 

ESCI310 

ESCI3S0 

ESC(3eC 


ESCI120 

ESCI13C 


ESCI170 

ESCIHOO 

ESCI19C 

ESCI5CC 

ESCI52C 

ESCI53? 


ESC(420 

ESCI430 
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nn 


TEMP»*l)«-TEMP»2)»EC0Mn7.N) 

TEMP»25)«TEHP«2)  •r.C0N«6,N) 

FM-.5*  <COS«P|*«  V2-V)/t V2-VI J )♦! .C) 

XV«l.O-ESTCON< I ,N)*V 
TI7>ESTC0NM7|N)-|  .333 
T2V-ESTC0NI 21 .N)/( I.C-X7)««2 

T3V«I  .0*XV*«  2.*T17*XV*ITI7*«  2«*FSTC0!>J(  I7  ,N)  >1  .6*47  )-ESf  CON  < 18.N)  ) ) 
TMV«T2V*T3V 

T‘tV>2.«IXV*ESTC0N(  la.tD-ESTCONI  17.N)*.3333) 

T5V«l  .2*XV*I  .E-l2*ESTC0NI6,fi)**2 
EC0N9V»EC0N«3tN)*<  CT2V*T‘»V*T3V)»*2)/T5V 
1 I VaESTCONI )4.N)*XV 

t2V>.5*l I.E-4*CSTC0N(4.h)*XVI**2/I 1.3-TIV) 

T3Val.*TlV*«l.*.5*TIV»«l  .-ESTCONCI 7,N)/ESTC0N« l6tN) ) )/3.0 
E0V«T2V*T3V*ESTC0tM23,N)*l  1 .♦ESTCONI  I 7 , fJ ) *X  V » ♦ESTCON « 22 , N ) 

IF  ITOLO  .LT.G.I  T0LD>3u3* 

TjaTOLD 

73  IF  ITJ  .r,T.  ECONIMtNI)  GO  TO  731 

VBaECONI ln.NI-Tj*(ECOIM  19.N)4EC0NI 2C  tNIkTJI 
TENP(7la>EC0N(  I 9 . tJ  I >2  . «ECON  I 20  • N ) • T J 
TEMP(8I*-2.*ECONI20,N) 

GO  TO  732 

731  VBaEC0N(6.H)/(Tj**.2S) 

TEMPI7)»-.25«EC0M4,N)/|  TJ**I.2S) 

TEMP(ei>S.*EC0N(6.Ni/( l6.«TJ*«2.25) 

732  ETA-VB/V 

TEMP(S)>(2.4>ETAI/( ( l.0-ETA)*«3) 

tv -TEMPI  I )«TJ-2.*CC0N|5.N)«TJ*Tj*TEPPCb) *TEHPI7) /W  ♦ E5TC0NI lOiN) 

. •|.E-l2*TEMPm 
XX«TJ/TMV 

EL*EOV*TEMPI IOI«TJ*.5»EcONI2.N)*Tj**2>ECONC3.N)»TMV*ECO!JlMtNI*.5» 

. EC0N9V^TEMP( I3I*TMV*I AL06I AUF*XX»I .OI-ECONC I5.N)-ALF*CXX-I.0) ) 
EJS«FM*EV4( I .-FMI*£L 
GO  TO  ICO 
C 

C •••  LIOUIO  REGION  ••• 

C 

79  |VFLAG(JI«32 

XX«T0LD/TEMP( 121 

TEMPI 23|aXX*l TEMPI Iw»*XA*TEMPI22) ) ♦TEMP  I I 3 ) • I ALOG I ALF»XX* I » I 
. -ALF*XXI 

EJSaTEMPI l2l*IEC0f)l3.HI*TEMPI |3)*l ALF-ECONI I5,N) )*TEMPI23) )♦ 

• TEMPI  1 1 l♦TEMPI6l 
GO  TO  ICC 

•••  VAPOR  REGION  ••• 

9C  IVFLAGIJ)-3S 

TEMPI9I«-TEMPI2I*CC0N| l/,NJ 
TEMPI25I«TEMPI2)*EC0N|4,N) 

TEMPI 17)«EC0MIS,N)*TEMP|25) 

TEMPI |8)a.5*TEMP||7l*TEMPI25) 

TTaTOLO 

91  IF  ITT  .GT.  ECONIII.N))  GO  TO  911 
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VB«EC0N(18,HI-TT*(EC0H( |».H»*EC0M«2C#N)»TT) 
TEMP{7)«-EC0N( lVtN)-2.*tC0N<25»MI*TT 
TEMP(e)a-2**tC0N(20tN) 

GO  TO  912 

911  V8»EC0N{ A,N)/(TT»».2SJ 
TEMP(7)«-.25«EC0N(6,NJ/(TT**1.2S) 
TEMP{a)«5.*EC0N(6.N)/( U.*TT**2.25) 

912  ETA-VB/V 

TEMP<5)«(2.S-ETA)/( < 1 .C-ETA»**3» 

EV-TEHPU  >*TT-2.*EC0N(5,NJ*TT*TT*TEMP(5) *TEHP(7l/V 

. *1 .E-I2*TEHP(«M 

EJS-EV 


* ESTCON(IOtN) 


10=  EJS«I>E>I2*EJS 
HETUKN 
END 


ESC22B 

ESC226 


APPENDIX  D 


SAMPLE  CALCULATION 

RIP  generation  and  restart  calculations  taken  from  the  Uranium 
study  are  presented  here  as  a sample  problem.  It  should  be  noted 
that  in  order  to  print  out  the  temperature  at  each  cycle,  the 
subroutine  MEDIT  as  listed  in  Reference  17  was  modified  in  the 
following  manner: 


The  line  MED  390 

PRINT  130,  J,  SMASS,  X(J),  RHO(J),  Q(J),  PT,  P(J),  E(J), 
C(J),  was  changed  to 

PRINT  130,  J,  SMASS,  X(J),  U(J),  RHO(J),  T(J),  PT,  P(J), 
E(J),  C(J), 

and  the  line  MED  680 

IHJ,  4X,  5HSMASS,  8X,  IHX,  lOX,  IHU,  lOX,  3HRH0,  8A,  IHQ, 
9X,  3HS+Q,  9X,  was  changed  to 

IHJ,  4X,  5HSMASS,  8X,  IHX,  lOX,  IHU,  lOX,  3HRH0,  8X, 

IHT,  9X,  3HS+Q,  9X, 
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